GCC Code Coverage Report

Directory: src/oiseau/
File: src/oiseau/dg/nodal/utils.cpp
Date: 2025-05-24 01:28:39
Exec Total Coverage
Lines: 524 524 100.0%
Functions: 38 38 100.0%
Branches: 592 1064 55.6%
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/utils.hpp"
8
9 #include <array>
10 #include <cmath>
11 #include <cstddef>
12 #include <numbers>
13 #include <utility>
14 #include <vector>
15 #include <xtensor-blas/xlinalg.hpp>
16 #include <xtensor/containers/xtensor.hpp>
17 #include <xtensor/core/xmath.hpp>
18 #include <xtensor/core/xoperation.hpp>
19 #include <xtensor/core/xtensor_forward.hpp>
20 #include <xtensor/generators/xbuilder.hpp>
21 #include <xtensor/misc/xmanipulation.hpp>
22 #include <xtensor/views/xslice.hpp>
23 #include <xtensor/views/xview.hpp>
24
25 #include "oiseau/utils/math.hpp"
26
27 namespace oiseau::dg::nodal::utils {
28
29 202 std::pair<xt::xarray<double>, xt::xarray<double>> jacobi_gq(unsigned n, double alpha, double beta) {
30
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 170 times.
 202 if (n == 0) {
31
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 xt::xarray<double> x = {(alpha - beta) / (alpha + beta + 2)};
32
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 xt::xarray<double> w = {2};
33 32 return make_pair(std::move(x), std::move(w));
34 32 }
35
1/2
✓ Branch 5 taken 170 times.
✗ Branch 6 not taken.
 170 xt::xarray<double> h1 = 2 * xt::linspace<double>(0, n, n + 1) + alpha + beta;
36
2/4
✓ Branch 1 taken 170 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 170 times.
✗ Branch 5 not taken.
 170 auto h1s = xt::view(h1, xt::range(0, n));
37 170 auto lins = xt::arange<double>(1, n + 1);
38 xt::xarray<double> m =
39
1/2
✓ Branch 6 taken 170 times.
✗ Branch 7 not taken.
 340 xt::diag(-1.0 / 2.0 * (std::pow(alpha, 2) - std::pow(beta, 2)) / (h1 + 2) / h1) +
40
1/2
✓ Branch 3 taken 170 times.
✗ Branch 4 not taken.
 510 xt::diag(2. / (h1s + 2) *
41 680 xt::sqrt(lins * (lins + alpha + beta) * (lins + alpha) * (lins + beta) /
42 850 (h1s + 1) / (h1s + 3)),
43
1/2
✓ Branch 2 taken 170 times.
✗ Branch 3 not taken.
 170 1);
44
3/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 158 times.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 170 if (alpha + beta == 0) m(0, 0) = 0;
45
1/2
✓ Branch 2 taken 170 times.
✗ Branch 3 not taken.
 170 m += xt::transpose(m);
46
1/2
✓ Branch 1 taken 170 times.
✗ Branch 2 not taken.
 170 auto [d, v] = xt::linalg::eigh(m);
47 170 double cte = std::pow(2, (alpha + beta + 1)) / (alpha + beta + 1) * std::tgamma(alpha + 1) *
48 170 std::tgamma(beta + 1) / std::tgamma(alpha + beta + 1);
49
3/6
✓ Branch 2 taken 170 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 170 times.
✗ Branch 8 not taken.
✓ Branch 13 taken 170 times.
✗ Branch 14 not taken.
 170 xt::xarray<double> w = cte * xt::square(xt::view(v, 0, xt::all())), x = d;
50
1/2
✓ Branch 1 taken 170 times.
✗ Branch 2 not taken.
 170 return {x, w};
51 170 }
52
53 196 xt::xarray<double> jacobi_gl(unsigned n, double alpha, double beta) {
54
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 177 times.
 196 if (n == 1) {
55
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 return xt::xarray<double>{-1.0, 1.0};
56 }
57
1/2
✓ Branch 1 taken 177 times.
✗ Branch 2 not taken.
 177 auto [x, _] = jacobi_gq(n - 2, alpha + 1.0, beta + 1.0);
58
2/4
✓ Branch 5 taken 177 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 177 times.
✗ Branch 9 not taken.
 177 return xt::concatenate(xt::xtuple(-xt::ones<double>({1}), x, xt::ones<double>({1})));
59 177 }
60
61 37 xt::xarray<double> vandermonde_1d(unsigned n, const xt::xarray<double> &r) {
62 37 std::array<size_t, 2> shape = {r.size(), n + 1};
63
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 37 xt::xtensor<double, 2> output(shape);
64
2/2
✓ Branch 0 taken 199 times.
✓ Branch 1 taken 37 times.
 236 for (size_t i = 0; i < n + 1; ++i) {
65
3/6
✓ Branch 1 taken 199 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 199 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 199 times.
✗ Branch 9 not taken.
 199 xt::view(output, xt::all(), i) = oiseau::utils::jacobi_p(i, 0.0, 0.0, r);
66 }
67
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 74 return output;
68 37 }
69
70 18 xt::xarray<double> grad_vandermonde_1d(unsigned n, const xt::xarray<double> &r) {
71 18 std::array<size_t, 2> shape = {r.size(), n + 1};
72
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 auto output = xt::xtensor<double, 2>::from_shape(shape);
73
2/2
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 18 times.
 138 for (size_t i = 0; i < n + 1; ++i) {
74
3/6
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 120 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 120 times.
✗ Branch 9 not taken.
 120 xt::view(output, xt::all(), i) = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, r);
75 }
76
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 return output;
77 18 }
78
79 17 xt::xarray<double> d_matrix_1d(const xt::xarray<double> &v, const xt::xarray<double> &gv) {
80
2/4
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 17 times.
✗ Branch 8 not taken.
 17 return xt::transpose(xt::linalg::solve(xt::transpose(v), xt::transpose(gv)));
81 }
82
83 9 xt::xarray<double> d_matrix_1d(unsigned n, const xt::xarray<double> &r) {
84
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto v = vandermonde_1d(n, r);
85
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto gv = grad_vandermonde_1d(n, r);
86
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 18 return d_matrix_1d(v, gv);
87 9 }
88
89 1 xt::xarray<double> generate_triangle_equidistant_nodes(unsigned n) {
90 1 std::size_t n_p = (n + 1) * (n + 2) / 2;
91
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 xt::xarray<double> l1 = xt::zeros<double>({n_p});
92
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 xt::xarray<double> l3 = xt::zeros<double>({n_p});
93
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 xt::xarray<double> out = xt::zeros<double>(xt::xarray<double>::shape_type{n_p, 3});
94
95 1 std::size_t idx = 0;
96
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 for (std::size_t i = 0; i < n + 1; ++i) {
97
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
 9 for (std::size_t j = 0; j < n + 1 - i; ++j, ++idx) {
98
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 l1[idx] = i * 1.0 / n;
99
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 l3[idx] = j * 1.0 / n;
100 }
101 }
102
103 #define l2 (1.0 - l1 - l3)
104
2/4
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
 1 xt::view(out, xt::all(), 0) = -l2 + l3;
105
2/4
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
 1 xt::view(out, xt::all(), 1) = (-l2 - l3 + 2.0 * l1) / std::numbers::sqrt3;
106 #undef l2
107 1 return out;
108 1 }
109
110 16 xt::xarray<double> warp_factor(unsigned n, const xt::xarray<double> &rout) {
111
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto lgl_r = oiseau::dg::nodal::utils::jacobi_gl(n, 0.0, 0.0);
112 16 auto r_eq = xt::linspace<double>(-1.0, 1.0, n + 1);
113
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 auto v_eq = oiseau::dg::nodal::utils::vandermonde_1d(n, r_eq);
114
2/4
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
 32 xt::xarray<double> p_mat = xt::zeros<double>(xt::xarray<double>::shape_type{n + 1, rout.size()});
115
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 16 times.
 80 for (unsigned i = 0; i < n + 1; ++i) {
116
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 xt::view(p_mat, i, xt::all()) = oiseau::utils::jacobi_p(i, 0.0, 0.0, rout);
117 }
118
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 auto l_mat = xt::linalg::solve(xt::transpose(v_eq), p_mat);
119
1/2
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 auto warp = xt::linalg::dot(xt::transpose(l_mat), (-r_eq + lgl_r));
120 16 auto zerof = xt::abs(rout) < 1.0 - 1e-10;
121 16 auto sf = 1.0 - xt::pow((zerof * rout), 2);
122
1/2
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 warp = warp / sf + warp * (zerof - 1);
123 16 return warp;
124 16 }
125
126 5 xt::xarray<double> conversion_equilateral_xy_to_rs(const xt::xarray<double> &coords) {
127
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto l_coords = xt::zeros_like(coords);
128
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto x = xt::col(coords, 0);
129
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto y = xt::col(coords, 1);
130 5 auto l1 = (std::numbers::sqrt3 * y + 1.0) / 3.0;
131 5 auto l2 = (-3.0 * x - std::numbers::sqrt3 * y + 2.0) / 6.0;
132 5 auto l3 = (3.0 * x - std::numbers::sqrt3 * y + 2.0) / 6.0;
133
2/4
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 xt::col(l_coords, 0) = -l2 + l3 - l1;
134
2/4
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 xt::col(l_coords, 1) = -l2 - (l3 - l1);
135 5 return l_coords;
136 5 }
137
138 67 xt::xarray<double> simplexp_2d(const xt::xarray<double> &ab, int i, int j) {
139
2/4
✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67 times.
✗ Branch 5 not taken.
 67 xt::xarray<double> a = xt::col(ab, 0);
140
2/4
✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67 times.
✗ Branch 5 not taken.
 67 xt::xarray<double> b = xt::col(ab, 1);
141
1/2
✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
 67 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
142
1/2
✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
 67 auto h2 = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
143
1/2
✓ Branch 6 taken 67 times.
✗ Branch 7 not taken.
 134 return std::numbers::sqrt2 * h1 * h2 * xt::pow(1 - b, i);
144 67 }
145
146 113 xt::xarray<double> tensorp_2d(const xt::xarray<double> &rs, int i, int j) {
147
2/4
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113 times.
✗ Branch 5 not taken.
 113 xt::xarray<double> r = xt::col(rs, 0);
148
2/4
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113 times.
✗ Branch 5 not taken.
 113 xt::xarray<double> s = xt::col(rs, 1);
149
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, r);
150
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 auto h2 = oiseau::utils::jacobi_p(j, 0.0, 0.0, s);
151
1/2
✓ Branch 2 taken 113 times.
✗ Branch 3 not taken.
 339 return h1 * h2;
152 113 }
153
154 216 xt::xarray<double> tensorp_3d(const xt::xarray<double> &rst, int i, int j, int k) {
155
2/4
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
 216 xt::xarray<double> r = xt::col(rst, 0);
156
2/4
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
 216 xt::xarray<double> s = xt::col(rst, 1);
157
2/4
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
 216 xt::xarray<double> t = xt::col(rst, 2);
158
159
1/2
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
 216 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, r); // P_i(r)
160
1/2
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
 216 auto h2 = oiseau::utils::jacobi_p(j, 0.0, 0.0, s); // P_j(s)
161
1/2
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
 216 auto h3 = oiseau::utils::jacobi_p(k, 0.0, 0.0, t); // P_k(t)
162
163
1/2
✓ Branch 3 taken 216 times.
✗ Branch 4 not taken.
 432 return h1 * h2 * h3;
164 216 }
165
166 36 xt::xarray<double> grad_tensorp_2d(const xt::xarray<double> &rs, int i, int j) {
167
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 xt::xarray<double> r = xt::col(rs, 0);
168
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 xt::xarray<double> s = xt::col(rs, 1);
169
170
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 auto pi_r = oiseau::utils::jacobi_p(i, 0.0, 0.0, r); // P_i(r)
171
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 auto dpi_r = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, r); // dP_i/dr
172
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 auto pj_s = oiseau::utils::jacobi_p(j, 0.0, 0.0, s); // P_j(s)
173
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 auto dpj_s = oiseau::utils::grad_jacobi_p(j, 0.0, 0.0, s); // dP_j/ds
174
175 36 auto dphidr = dpi_r * pj_s; // ∂φ/∂r
176 36 auto dphids = pi_r * dpj_s; // ∂φ/∂s
177
178
2/4
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
 72 return xt::stack(xt::xtuple(dphidr, dphids), 1); // shape: [N, 2]
179 36 }
180
181 125 xt::xarray<double> grad_tensorp_3d(const xt::xarray<double> &rst, int i, int j, int k) {
182
2/4
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
 125 xt::xarray<double> r = xt::col(rst, 0);
183
2/4
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
 125 xt::xarray<double> s = xt::col(rst, 1);
184
2/4
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
 125 xt::xarray<double> t = xt::col(rst, 2);
185
186
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto Pi_r = oiseau::utils::jacobi_p(i, 0.0, 0.0, r); // P_i(r)
187
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto Pj_s = oiseau::utils::jacobi_p(j, 0.0, 0.0, s); // P_j(s)
188
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto Pk_t = oiseau::utils::jacobi_p(k, 0.0, 0.0, t); // P_k(t)
189
190
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto dPi_r = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, r); // dP_i/dr
191
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto dPj_s = oiseau::utils::grad_jacobi_p(j, 0.0, 0.0, s); // dP_j/ds
192
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto dPk_t = oiseau::utils::grad_jacobi_p(k, 0.0, 0.0, t); // dP_k/dt
193
194 125 auto dphidr = dPi_r * Pj_s * Pk_t;
195 125 auto dphids = Pi_r * dPj_s * Pk_t;
196 125 auto dphidt = Pi_r * Pj_s * dPk_t;
197
198
2/4
✓ Branch 2 taken 125 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 125 times.
✗ Branch 6 not taken.
 250 return xt::stack(xt::xtuple(dphidr, dphids, dphidt), 1); // shape: [N, 3]
199 125 }
200
201 21 xt::xarray<double> grad_simplexp_2d(const xt::xarray<double> &ab, int i, int j) {
202
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 21 xt::xarray<double> a = xt::col(ab, 0);
203
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 21 xt::xarray<double> b = xt::col(ab, 1);
204
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> fa = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
205
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> dfa = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, a);
206
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> gb = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
207
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> dgb = oiseau::utils::grad_jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
208
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 xt::xarray<double> dmodedr = dfa * gb;
209
1/2
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 21 xt::xarray<double> dmodeds = dfa * (gb * (0.5 * (1 + a)));
210
1/2
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 21 xt::xarray<double> tmp = dgb * xt::pow(0.5 * (1 - b), i);
211
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 8 times.
 21 if (i > 0) {
212
1/2
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
 13 dmodedr *= xt::pow(0.5 * (1 - b), i - 1);
213
1/2
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
 13 dmodeds *= xt::pow(0.5 * (1 - b), i - 1);
214
1/2
✓ Branch 6 taken 13 times.
✗ Branch 7 not taken.
 13 tmp += -0.5 * i * gb * xt::pow(0.5 * (1 - b), i - 1);
215 }
216
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 dmodeds = dmodeds + fa * tmp;
217
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 dmodedr = pow(2.0, i + 0.5) * dmodedr;
218
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 dmodeds = pow(2.0, i + 0.5) * dmodeds;
219
2/4
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 42 return xt::stack(xt::xtuple(dmodedr, dmodeds), 1);
220 21 }
221
222 5 xt::xarray<double> vandermonde_2d(unsigned n, const xt::xarray<double> &rs) {
223
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto ab = conversion_rs_to_ab(rs);
224 5 std::array<size_t, 2> shape = {rs.shape()[0], ((n + 1) * (n + 2)) / 2};
225
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 xt::xtensor<double, 2> output(shape);
226 5 int index = 0;
227
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 5 times.
 28 for (unsigned i = 0; i <= n; ++i) {
228
2/2
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 23 times.
 89 for (unsigned j = 0; j <= n - i; ++j, ++index) {
229
3/6
✓ Branch 1 taken 66 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 66 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 66 times.
✗ Branch 8 not taken.
 66 xt::col(output, index) = simplexp_2d(ab, i, j);
230 }
231 }
232
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 return output;
233 5 }
234
235 4 xt::xarray<double> vandermonde_2d_tensor(unsigned n, const xt::xarray<double> &rs) {
236 4 std::array<size_t, 2> shape = {rs.shape()[0], (n + 1) * (n + 1)};
237
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 xt::xtensor<double, 2> output(shape);
238 4 int index = 0;
239
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 4 times.
 25 for (unsigned i = 0; i <= n; ++i) {
240
2/2
✓ Branch 0 taken 113 times.
✓ Branch 1 taken 21 times.
 134 for (unsigned j = 0; j <= n; ++j, ++index) {
241
5/10
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 113 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 113 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 113 times.
✗ Branch 14 not taken.
 113 xt::col(output, index) = xt::col(output, index) = tensorp_2d(rs, i, j);
242 }
243 }
244
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 return output;
245 4 }
246
247 2 xt::xarray<double> grad_vandermonde_2d(unsigned n, const xt::xarray<double> &rs) {
248
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto ab = conversion_rs_to_ab(rs);
249 2 std::array<size_t, 3> shape = {rs.shape()[0], ((n + 1) * (n + 2)) / 2, 2};
250
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 xt::xtensor<double, 3> output(shape);
251 2 int index = 0;
252
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for (unsigned i = 0; i <= n; ++i) {
253
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 8 times.
 28 for (unsigned j = 0; j <= n - i; ++j, ++index) {
254
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 auto tmp = grad_simplexp_2d(ab, i, j);
255
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
 20 xt::view(output, xt::all(), index, 0) = xt::col(tmp, 0);
256
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
 20 xt::view(output, xt::all(), index, 1) = xt::col(tmp, 1);
257 20 }
258 }
259
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 return output;
260 2 }
261
262 1 xt::xarray<double> grad_vandermonde_2d_tensor(unsigned n, const xt::xarray<double> &rs) {
263 1 std::array<size_t, 3> shape = {rs.shape()[0], (n + 1) * (n + 1), 2};
264
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 xt::xtensor<double, 3> output(shape);
265 1 int index = 0;
266
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 times.
 7 for (unsigned i = 0; i <= n; ++i) {
267
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 6 times.
 42 for (unsigned j = 0; j <= n; ++j, ++index) {
268
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 auto tmp = grad_tensorp_2d(rs, i, j);
269
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 36 times.
✗ Branch 9 not taken.
 36 xt::view(output, xt::all(), index, 0) = xt::col(tmp, 0);
270
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 36 times.
✗ Branch 9 not taken.
 36 xt::view(output, xt::all(), index, 1) = xt::col(tmp, 1);
271 36 }
272 }
273
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 return output;
274 1 }
275
276 2 xt::xarray<double> d_matrix_2d(const xt::xarray<double> &v, const xt::xarray<double> &gv) {
277
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
278
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto gvs = xt::view(gv, xt::all(), xt::all(), 1);
279
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto dr = xt::transpose(xt::linalg::solve(xt::transpose(v), xt::transpose(gvr)));
280
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto ds = xt::transpose(xt::linalg::solve(xt::transpose(v), xt::transpose(gvs)));
281
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 4 return xt::stack(xt::xtuple(dr, ds), 2);
282 2 }
283
284 9 xt::xarray<double> conversion_rs_to_ab(const xt::xarray<double> &rs) {
285
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto ab = xt::zeros_like(rs);
286
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto r = xt::col(rs, 0);
287
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto s = xt::col(rs, 1);
288
2/2
✓ Branch 2 taken 145 times.
✓ Branch 3 taken 9 times.
 154 for (unsigned i = 0; i < rs.shape()[0]; ++i) {
289
6/10
✓ Branch 1 taken 145 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 136 times.
✓ Branch 4 taken 9 times.
✓ Branch 6 taken 136 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 136 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 145 times.
✗ Branch 13 not taken.
 145 ab(i, 0) = (s[i] == 1) ? -1 : 2 * (1 + r[i]) / (1 - s[i]) - 1;
290 }
291
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 xt::col(ab, 1) = s;
292 9 return ab;
293 9 }
294
295 5 xt::xarray<double> generate_triangle_nodes(unsigned n) {
296 5 constexpr auto alp_opt =
297 std::to_array({0.0000, 0.0000, std::numbers::sqrt2, 0.1001, 0.2751, 0.9800, 1.0999, 1.2832,
298 1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258});
299
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 double alpha = (n < 16) ? alp_opt[n - 1] : 5.0 / 3.0;
300 5 std::size_t n_p = (n + 1) * (n + 2) / 2;
301
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 xt::xarray<double> l1 = xt::zeros<double>({n_p});
302
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 xt::xarray<double> l3 = xt::zeros<double>({n_p});
303
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
 10 xt::xarray<double> out = xt::zeros<double>(xt::xarray<double>::shape_type{n_p, 2});
304 5 std::size_t idx = 0;
305
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
 25 for (std::size_t i = 0; i < n + 1; ++i) {
306
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 20 times.
 70 for (std::size_t j = 0; j < n + 1 - i; ++j, ++idx) {
307
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 l1[idx] = i * 1.0 / n;
308
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 l3[idx] = j * 1.0 / n;
309 }
310 }
311
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 xt::xarray<double> l2 = (1.0 - l1 - l3);
312 5 auto blend1 = 4 * l2 * l3;
313 5 auto blend2 = 4 * l1 * l3;
314 5 auto blend3 = 4 * l1 * l2;
315
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 xt::xarray<double> e1 = -l2 + l3;
316
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 xt::xarray<double> e2 = -l3 + l1;
317
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 xt::xarray<double> e3 = -l1 + l2;
318
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto warpf1 = oiseau::dg::nodal::utils::warp_factor(n, e1);
319
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto warpf2 = oiseau::dg::nodal::utils::warp_factor(n, e2);
320
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto warpf3 = oiseau::dg::nodal::utils::warp_factor(n, e3);
321 5 auto warp1 = blend1 * warpf1 * xt::square(1 + (alpha * l1));
322 5 auto warp2 = blend2 * warpf2 * xt::square(1 + (alpha * l2));
323 5 auto warp3 = blend3 * warpf3 * xt::square(1 + (alpha * l3));
324
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto &&x = xt::col(out, 0);
325
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto &&y = xt::col(out, 1);
326 5 auto pi = xt::numeric_constants<double>::PI;
327
2/4
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
 5 xt::col(out, 0) = -l2 + l3;
328
2/4
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
 5 xt::col(out, 1) = (-l2 - l3 + 2.0 * l1) / std::numbers::sqrt3;
329
1/2
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 x = x + 1 * warp1 + std::cos(2 * pi / 3) * warp2 + std::cos(4 * pi / 3) * warp3;
330
1/2
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 y = y + 0 * warp1 + std::sin(2 * pi / 3) * warp2 + std::sin(4 * pi / 3) * warp3;
331 5 return out;
332 5 }
333
334 1 xt::xarray<double> generate_quadrilateral_nodes(unsigned order) {
335
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 xt::xarray<double> r1d = jacobi_gl(order, 0.0, 0.0);
336 1 unsigned Np1 = order + 1;
337 1 unsigned Np_quad = Np1 * Np1;
338
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 std::vector<std::size_t> shape = {static_cast<std::size_t>(Np_quad), 2ul};
339
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 xt::xarray<double> rs_nodes = xt::zeros<double>(shape);
340
341 1 unsigned k = 0;
342
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 times.
 7 for (unsigned i = 0; i < Np1; ++i) {
343
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 6 times.
 42 for (unsigned j = 0; j < Np1; ++j) {
344
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 rs_nodes(k, 0) = r1d(j);
345
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 rs_nodes(k, 1) = r1d(i);
346 36 k++;
347 }
348 }
349 1 return rs_nodes;
350 1 }
351
352 1 xt::xarray<double> generate_hexahedron_nodes(unsigned order) {
353
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 xt::xarray<double> r1d = jacobi_gl(order, 0.0, 0.0);
354
355 1 unsigned Np1 = order + 1;
356 1 std::size_t Np_hex = Np1 * Np1 * Np1;
357
358
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 std::vector<std::size_t> shape = {Np_hex, 3ul};
359
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 xt::xarray<double> rst_nodes = xt::zeros<double>(shape);
360
361 1 unsigned k_idx = 0;
362
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 times.
 6 for (unsigned k_loop = 0; k_loop < Np1; ++k_loop) {
363
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 5 times.
 30 for (unsigned i_loop = 0; i_loop < Np1; ++i_loop) {
364
2/2
✓ Branch 0 taken 125 times.
✓ Branch 1 taken 25 times.
 150 for (unsigned j_loop = 0; j_loop < Np1; ++j_loop) {
365
2/4
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
 125 rst_nodes(k_idx, 0) = r1d(j_loop);
366
2/4
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
 125 rst_nodes(k_idx, 1) = r1d(i_loop);
367
2/4
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
 125 rst_nodes(k_idx, 2) = r1d(k_loop);
368 125 k_idx++;
369 }
370 }
371 }
372 1 return rst_nodes;
373 1 }
374
375 7 xt::xarray<double> generate_tetrahedron_equidistant_nodes(unsigned n) {
376 7 auto n_p = ((n + 1) * (n + 2) * (n + 3)) / 6;
377
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 auto shape = xt::xarray<double>::shape_type{n_p, 3};
378
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
 7 xt::xarray<double> output = xt::zeros<double>(shape);
379 7 std::size_t sk = 0;
380
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 7 times.
 31 for (std::size_t i = 0; i < n + 1; ++i) {
381
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 24 times.
 78 for (std::size_t j = 0; j < n + 1 - i; ++j) {
382
2/2
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 54 times.
 154 for (std::size_t k = 0; k < n + 1 - i - j; ++k, ++sk) {
383
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 output(sk, 0) = -1 + k * 2.0 / n;
384
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 output(sk, 1) = -1 + j * 2.0 / n;
385
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 output(sk, 2) = -1 + i * 2.0 / n;
386 }
387 }
388 }
389 7 return output;
390 7 }
391
392 76 xt::xarray<double> eval_warp(int p, const xt::xarray<double> &xnodes,
393 const xt::xarray<double> &xout) {
394
1/2
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
 76 xt::xarray<double> warp = xt::zeros_like(xout);
395
1/2
✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
 76 xt::xarray<double> xeq = xt::zeros<double>({p + 1});
396
3/4
✓ Branch 1 taken 268 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 268 times.
✓ Branch 4 taken 76 times.
 344 for (int i = 0; i <= p; i++) xeq(i) = -1.0 + 2.0 * (p - i) / p;
397
2/2
✓ Branch 0 taken 268 times.
✓ Branch 1 taken 76 times.
 344 for (int i = 0; i <= p; i++) {
398
4/8
✓ Branch 1 taken 268 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 268 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 268 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 268 times.
✗ Branch 12 not taken.
 268 xt::xarray<double> d = xt::ones_like(xout) * (xnodes(i) - xeq(i));
399
2/2
✓ Branch 0 taken 428 times.
✓ Branch 1 taken 268 times.
 696 for (int j = 1; j < p; j++) {
400
6/10
✓ Branch 0 taken 312 times.
✓ Branch 1 taken 116 times.
✓ Branch 3 taken 312 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 312 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 312 times.
✗ Branch 10 not taken.
✓ Branch 15 taken 312 times.
✗ Branch 16 not taken.
 428 if (i != j) d = d * (xout - xeq(j)) / (xeq(i) - xeq(j));
401 }
402
5/8
✓ Branch 0 taken 192 times.
✓ Branch 1 taken 76 times.
✓ Branch 3 taken 192 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 192 times.
✗ Branch 7 not taken.
✓ Branch 11 taken 192 times.
✗ Branch 12 not taken.
 268 if (i != 0) d = -d / (xeq(i) - xeq(0));
403
5/8
✓ Branch 0 taken 192 times.
✓ Branch 1 taken 76 times.
✓ Branch 3 taken 192 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 192 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 192 times.
✗ Branch 11 not taken.
 268 if (i != p) d = d / (xeq(i) - xeq(p));
404
1/2
✓ Branch 1 taken 268 times.
✗ Branch 2 not taken.
 268 warp += d;
405 268 }
406 76 return warp;
407 76 }
408
409 25 xt::xarray<double> eval_shift(int p, double pval, const xt::xarray<double> &l1,
410 const xt::xarray<double> &l2, const xt::xarray<double> &l3) {
411
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 auto gaussX = -jacobi_gl(p, 0.0, 0.0);
412 25 auto blend1 = l2 * l3;
413 25 auto blend2 = l1 * l3;
414 25 auto blend3 = l1 * l2;
415
3/6
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 25 times.
✗ Branch 10 not taken.
 25 auto warpfactor1 = 4 * eval_warp(p, gaussX, l3 - 1 * l2);
416
3/6
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 25 times.
✗ Branch 10 not taken.
 25 auto warpfactor2 = 4 * eval_warp(p, gaussX, l1 - 1 * l3);
417
3/6
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 25 times.
✗ Branch 10 not taken.
 25 auto warpfactor3 = 4 * eval_warp(p, gaussX, l2 - 1 * l1);
418 25 auto warp1 = blend1 * warpfactor1 * (1 + xt::pow(pval * l1, 2));
419 25 auto warp2 = blend2 * warpfactor2 * (1 + xt::pow(pval * l2, 2));
420 25 auto warp3 = blend3 * warpfactor3 * (1 + xt::pow(pval * l3, 2));
421 25 auto pi = xt::numeric_constants<double>::PI;
422 25 auto dx = 1 * warp1 + cos(2.0 * pi / 3.0) * warp2 + cos(4.0 * pi / 3.0) * warp3;
423 25 auto dy = 0 * warp1 + sin(2.0 * pi / 3.0) * warp2 + sin(4.0 * pi / 3.0) * warp3;
424
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 50 return xt::stack(xt::xtuple(dx, dy), 1);
425 25 }
426
427 24 xt::xarray<double> warp_shift_face_3d(int p, double pval, const xt::xarray<double> &l2,
428 const xt::xarray<double> &l3, const xt::xarray<double> &l4) {
429 24 return eval_shift(p, pval, l2, l3, l4);
430 }
431
432 6 xt::xarray<double> generate_tetrahedron_nodes(unsigned p) {
433 xt::xarray<double> alphastore = {
434 0, 0, 0, 0.1002, 1.1332, 1.5608, 1.3413, 1.2577,
435 1.1603, 1.10153, 0.6080, 0.4523, 0.8856, 0.8717, 0.9655,
436
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 };
437
2/4
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
 6 double alpha = p > 14 ? 1 : alphastore[p];
438 6 double tol = 1e-10;
439
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 auto rst = generate_tetrahedron_equidistant_nodes(p);
440
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 6 auto r = xt::col(rst, 0), s = xt::col(rst, 1), t = xt::col(rst, 2);
441
3/6
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 20 taken 6 times.
✗ Branch 21 not taken.
 6 xt::xarray<double> l1 = (1 + t) / 2, l2 = (1 + s) / 2, l3 = -(1 + r + s + t) / 2,
442
1/2
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
 6 l4 = (1 + r) / 2;
443
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 xt::xtensor<double, 1> v1 = {-1.0, -1.0 / std::numbers::sqrt3, -1.0 / std::sqrt(6)};
444
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 xt::xtensor<double, 1> v2 = {1.0, -1.0 / std::numbers::sqrt3, -1.0 / std::sqrt(6)};
445
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 xt::xtensor<double, 1> v3 = {0.0, 2.0 / std::numbers::sqrt3, -1.0 / std::sqrt(6)};
446
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 xt::xtensor<double, 1> v4 = {0.0, 0.0, 3.0 / std::sqrt(6)};
447
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 xt::xtensor<double, 2> t1 = xt::zeros<double>({4, 3});
448
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 xt::xtensor<double, 2> t2 = xt::zeros<double>({4, 3});
449
2/4
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
 6 xt::row(t1, 0) = v2 - 1 * v1;
450
2/4
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
 6 xt::row(t1, 1) = v2 - 1 * v1;
451
2/4
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
 6 xt::row(t1, 2) = v3 - 1 * v2;
452
2/4
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
 6 xt::row(t1, 3) = v3 - 1 * v1;
453
2/4
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 6 xt::row(t2, 0) = v3 - 0.5 * (v1 + v2);
454
2/4
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 6 xt::row(t2, 1) = v4 - 0.5 * (v1 + v2);
455
2/4
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 6 xt::row(t2, 2) = v4 - 0.5 * (v2 + v3);
456
2/4
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 6 xt::row(t2, 3) = v4 - 0.5 * (v1 + v3);
457
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 6 times.
 30 for (int i = 0; i < 4; ++i) {
458
4/8
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
 24 xt::row(t1, i) /= xt::linalg::norm(xt::row(t1, i));
459
4/8
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
 24 xt::row(t2, i) /= xt::linalg::norm(xt::row(t2, i));
460 }
461
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 12 auto xyz = xt::linalg::outer(l3, v1) + xt::linalg::outer(l4, v2) + xt::linalg::outer(l2, v3) +
462
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 18 xt::linalg::outer(l1, v4);
463
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 6 xt::xarray<double> shift = xt::zeros<double>(xyz.shape());
464
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 6 times.
 30 for (int face = 0; face < 4; ++face) {
465
4/8
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
 24 xt::xarray<double> la, lb, lc, ld;
466
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 18 times.
 24 if (face == 0) {
467
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 la = l1;
468
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lb = l2;
469
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lc = l3;
470
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 ld = l4;
471
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 12 times.
 18 } else if (face == 1) {
472
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 la = l2;
473
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lb = l1;
474
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lc = l3;
475
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 ld = l4;
476
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 6 times.
 12 } else if (face == 2) {
477
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 la = l3;
478
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lb = l1;
479
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lc = l4;
480
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 ld = l2;
481
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 } else if (face == 3) {
482
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 la = l4;
483
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lb = l1;
484
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 lc = l3;
485
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 ld = l2;
486 }
487
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 auto warp = warp_shift_face_3d(p, alpha, lb, lc, ld);
488
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 auto warp1 = xt::col(warp, 0);
489
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 auto warp2 = xt::col(warp, 1);
490
1/2
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
 24 xt::xarray<double> blend = lb * lc * ld;
491
1/2
✓ Branch 9 taken 24 times.
✗ Branch 10 not taken.
 24 xt::xarray<double> denom = (lb + 0.5 * la) * (lc + 0.5 * la) * (ld + 0.5 * la);
492
2/4
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
 24 auto ids = xt::flatten_indices(xt::argwhere(denom > tol));
493
4/8
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 24 times.
✗ Branch 14 not taken.
 96 xt::view(blend, xt::keep(ids)) = (1 + xt::pow(alpha * xt::view(la, xt::keep(ids)), 2)) *
494
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 96 xt::view(blend, xt::keep(ids)) /
495
3/6
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 24 times.
✗ Branch 9 not taken.
 96 xt::view(denom, xt::keep(ids));
496
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
 48 shift = shift + xt::linalg::outer(blend * warp1, xt::row(t1, face)) +
497
3/6
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 24 times.
✗ Branch 10 not taken.
 72 xt::linalg::outer(blend * warp2, xt::row(t2, face));
498 24 auto mask = (la < tol) && (((lb > tol) + (lc > tol) + (ld > tol)) < 3);
499
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 24 ids = xt::flatten_indices(xt::argwhere(mask));
500
2/4
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
 48 xt::view(shift, xt::keep(ids), xt::all()) =
501
4/8
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
 48 xt::linalg::outer(xt::view(warp1, xt::keep(ids)), xt::row(t1, face)) +
502
5/10
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 24 times.
✗ Branch 15 not taken.
 96 xt::linalg::outer(xt::view(warp2, xt::keep(ids)), xt::row(t2, face));
503 24 }
504 6 xyz = xyz + shift;
505
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 return xyz;
506 6 }
507
508 5 xt::xarray<double> conversion_equilateral_xyz_to_rst(const xt::xarray<double> &coords) {
509 5 const double sqrt3 = std::numbers::sqrt3;
510 5 const double sqrt6 = std::sqrt(6.0);
511
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 xt::xtensor<double, 1> v1 = {-1, -1 / sqrt3, -1 / sqrt6};
512
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 xt::xtensor<double, 1> v2 = {1, -1 / sqrt3, -1 / sqrt6};
513
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 xt::xtensor<double, 1> v3 = {0, 2 / sqrt3, -1 / sqrt6};
514
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 xt::xtensor<double, 1> v4 = {0, 0 / sqrt3, 3 / sqrt6};
515 5 auto rhs = xt::transpose(coords - 0.5 * (v2 + v3 + v4 - v1));
516
1/2
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
 5 auto lhs = 0.5 * xt::stack(xt::xtuple(v2 - v1, v3 - v1, v4 - v1), 1);
517
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto rst = xt::linalg::solve(lhs, rhs);
518
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 10 return xt::transpose(rst);
519 5 }
520
521 142 xt::xarray<double> simplexp_3d(const xt::xarray<double> &abc, int i, int j, int k) {
522
2/4
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 142 times.
✗ Branch 5 not taken.
 142 xt::xarray<double> a = xt::col(abc, 0);
523
2/4
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 142 times.
✗ Branch 5 not taken.
 142 xt::xarray<double> b = xt::col(abc, 1);
524
2/4
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 142 times.
✗ Branch 5 not taken.
 142 xt::xarray<double> c = xt::col(abc, 2);
525
1/2
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
 142 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
526
1/2
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
 142 auto h2 = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
527
1/2
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
 142 auto h3 = oiseau::utils::jacobi_p(k, 2.0 * (i + j) + 2.0, 0.0, c);
528
1/2
✓ Branch 10 taken 142 times.
✗ Branch 11 not taken.
 284 return 2 * std::numbers::sqrt2 * h1 * h2 * xt::pow(1 - b, i) * h3 * xt::pow(1 - c, i + j);
529 142 }
530
531 9 xt::xarray<double> conversion_rst_to_abc(const xt::xarray<double> &rst) {
532
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto abc = xt::zeros_like(rst);
533
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto &&r = xt::col(rst, 0);
534
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto &&s = xt::col(rst, 1);
535
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto &&t = xt::col(rst, 2);
536
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto &&a = xt::col(abc, 0);
537
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto &&b = xt::col(abc, 1);
538
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 auto &&c = xt::col(abc, 2);
539
2/2
✓ Branch 2 taken 350 times.
✓ Branch 3 taken 9 times.
 359 for (unsigned n = 0; n < rst.shape()[0]; ++n) {
540
4/6
✓ Branch 1 taken 350 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 350 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 321 times.
✓ Branch 7 taken 29 times.
 350 if ((s[n] + t[n]) != 0) {
541
4/8
✓ Branch 1 taken 321 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 321 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 321 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 321 times.
✗ Branch 11 not taken.
 321 a[n] = 2 * ((1 + r[n]) / (-s[n] - t[n])) - 1;
542 } else {
543
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 29 a[n] = -1;
544 }
545
3/4
✓ Branch 1 taken 350 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 341 times.
✓ Branch 4 taken 9 times.
 350 if (t[n] != 1) {
546
3/6
✓ Branch 1 taken 341 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 341 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 341 times.
✗ Branch 8 not taken.
 341 b[n] = 2 * (1 + s[n]) / (1 - t[n]) - 1;
547 } else {
548
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 b[n] = -1;
549 }
550 }
551
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 c = t;
552 9 return abc;
553 9 }
554
555 5 xt::xarray<double> vandermonde_3d(unsigned n, const xt::xarray<double> &rst) {
556
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 auto abc = conversion_rst_to_abc(rst);
557 5 std::array<size_t, 2> shape = {rst.shape()[0], ((n + 1) * (n + 2) * (n + 3)) / 6};
558
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 xt::xtensor<double, 2> output(shape);
559 5 int index = 0;
560
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 5 times.
 27 for (unsigned i = 0; i <= n; ++i) {
561
2/2
✓ Branch 0 taken 62 times.
✓ Branch 1 taken 22 times.
 84 for (unsigned j = 0; j <= n - i; ++j) {
562
2/2
✓ Branch 0 taken 141 times.
✓ Branch 1 taken 62 times.
 203 for (unsigned k = 0; k <= n - i - j; ++k, ++index) {
563
3/6
✓ Branch 1 taken 141 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 141 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 141 times.
✗ Branch 8 not taken.
 141 xt::col(output, index) = simplexp_3d(abc, i, j, k);
564 }
565 }
566 }
567
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 return output;
568 5 }
569
570 3 xt::xarray<double> vandermonde_3d_tensor(unsigned n, const xt::xarray<double> &rst) {
571 3 std::array<size_t, 2> shape = {rst.shape()[0], (n + 1) * (n + 1) * (n + 1)};
572
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 xt::xtensor<double, 2> output(shape);
573 3 int index = 0;
574
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 3 times.
 15 for (unsigned i = 0; i <= n; ++i) {
575
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 12 times.
 62 for (unsigned j = 0; j <= n; ++j) {
576
2/2
✓ Branch 0 taken 216 times.
✓ Branch 1 taken 50 times.
 266 for (unsigned k = 0; k <= n; ++k, ++index) {
577
3/6
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 216 times.
✗ Branch 8 not taken.
 216 xt::col(output, index) = tensorp_3d(rst, i, j, k);
578 }
579 }
580 }
581
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 return output;
582 3 }
583
584 21 xt::xarray<double> grad_simplexp_3d(const xt::xarray<double> &abc, int i, int j, int k) {
585
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 21 xt::xarray<double> a = xt::col(abc, 0);
586
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 21 xt::xarray<double> b = xt::col(abc, 1);
587
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 21 xt::xarray<double> c = xt::col(abc, 2);
588
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> fa = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
589
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> dfa = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, a);
590
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> gb = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
591
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> dgb = oiseau::utils::grad_jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
592
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> hc = oiseau::utils::jacobi_p(k, 2.0 * (i + j) + 2.0, 0.0, c);
593
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 xt::xarray<double> dhc = oiseau::utils::grad_jacobi_p(k, 2.0 * (i + j) + 2.0, 0.0, c);
594
595
4/8
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 21 times.
✗ Branch 11 not taken.
 21 xt::xarray<double> v3dr, v3ds, v3dt, tmp;
596
597
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 v3dr = dfa * (gb * hc);
598
3/4
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 12 times.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
 21 if (i > 0) v3dr = v3dr * xt::pow(0.5 * (1 - b), i - 1);
599
3/4
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 6 times.
✓ Branch 7 taken 15 times.
✗ Branch 8 not taken.
 21 if (i + j > 0) v3dr = v3dr * xt::pow(0.5 * (1 - c), i + j - 1);
600
601
1/2
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 21 v3ds = 0.5 * (1 + a) * v3dr;
602
1/2
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 21 tmp = dgb * xt::pow(0.5 * (1 - b), i);
603
3/4
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 12 times.
✓ Branch 9 taken 9 times.
✗ Branch 10 not taken.
 21 if (i > 0) tmp = tmp + (-0.5 * i) * gb * xt::pow(0.5 * (1 - b), i - 1);
604
3/4
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 6 times.
✓ Branch 7 taken 15 times.
✗ Branch 8 not taken.
 21 if (i + j > 0) tmp = tmp * xt::pow(0.5 * (1 - c), i + j - 1);
605
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 tmp = fa * (tmp * hc);
606
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 v3ds = v3ds + tmp;
607
608
1/2
✓ Branch 8 taken 21 times.
✗ Branch 9 not taken.
 21 v3dt = 0.5 * (1 + a) * v3dr + 0.5 * (1 + b) * tmp;
609
1/2
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 21 tmp = dhc * xt::pow(0.5 * (1 - c), i + j);
610
3/4
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 6 times.
✓ Branch 9 taken 15 times.
✗ Branch 10 not taken.
 21 if (i + j > 0) tmp = tmp - 0.5 * (i + j) * hc * xt::pow(0.5 * (1 - c), (i + j - 1));
611
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 tmp = fa * gb * tmp;
612
1/2
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 21 tmp = tmp * xt::pow(0.5 * (1 - b), i);
613
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 v3dt = v3dt + tmp;
614
615
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 v3dr = v3dr * std::pow(2, 2 * i + j + 1.5);
616
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 v3ds = v3ds * std::pow(2, 2 * i + j + 1.5);
617
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 21 v3dt = v3dt * std::pow(2, 2 * i + j + 1.5);
618
619
2/4
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 42 return xt::stack(xt::xtuple(v3dr, v3ds, v3dt), 1);
620 21 }
621
622 2 xt::xarray<double> grad_vandermonde_3d(unsigned n, const xt::xarray<double> &rst) {
623
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto abc = conversion_rst_to_abc(rst);
624 2 std::array<size_t, 3> shape = {rst.shape()[0], ((n + 1) * (n + 2) * (n + 3)) / 6, 3};
625
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 xt::xtensor<double, 3> output(shape);
626 2 int index = 0;
627
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 for (unsigned i = 0; i <= n; ++i) {
628
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 6 times.
 18 for (unsigned j = 0; j <= n - i; ++j) {
629
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 12 times.
 32 for (unsigned k = 0; k <= n - i - j; ++k, ++index) {
630
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 auto tmp = grad_simplexp_3d(abc, i, j, k);
631
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
 20 xt::view(output, xt::all(), index, 0) = xt::col(tmp, 0);
632
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
 20 xt::view(output, xt::all(), index, 1) = xt::col(tmp, 1);
633
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
 20 xt::view(output, xt::all(), index, 2) = xt::col(tmp, 2);
634 20 }
635 }
636 }
637
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 return output;
638 2 }
639
640 1 xt::xarray<double> grad_vandermonde_3d_tensor(unsigned n, const xt::xarray<double> &rst) {
641 1 std::array<size_t, 3> shape = {rst.shape()[0], (n + 1) * (n + 1) * (n + 1), 3};
642
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 xt::xtensor<double, 3> output(shape);
643 1 int index = 0;
644
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 times.
 6 for (unsigned i = 0; i <= n; ++i) {
645
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 5 times.
 30 for (unsigned j = 0; j <= n; ++j) {
646
2/2
✓ Branch 0 taken 125 times.
✓ Branch 1 taken 25 times.
 150 for (unsigned k = 0; k <= n; ++k, ++index) {
647
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 auto tmp = grad_tensorp_3d(rst, i, j, k);
648
3/6
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 125 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 125 times.
✗ Branch 9 not taken.
 125 xt::view(output, xt::all(), index, 0) = xt::col(tmp, 0);
649
3/6
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 125 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 125 times.
✗ Branch 9 not taken.
 125 xt::view(output, xt::all(), index, 1) = xt::col(tmp, 1);
650
3/6
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 125 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 125 times.
✗ Branch 9 not taken.
 125 xt::view(output, xt::all(), index, 2) = xt::col(tmp, 2);
651 125 }
652 }
653 }
654
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 return output;
655 1 }
656
657 2 xt::xarray<double> d_matrix_3d(const xt::xarray<double> &v, const xt::xarray<double> &gv) {
658
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
659
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto gvs = xt::view(gv, xt::all(), xt::all(), 1);
660
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto gvt = xt::view(gv, xt::all(), xt::all(), 2);
661
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto dr = xt::transpose(xt::linalg::solve(xt::transpose(v), xt::transpose(gvr)));
662
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto ds = xt::transpose(xt::linalg::solve(xt::transpose(v), xt::transpose(gvs)));
663
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 auto dt = xt::transpose(xt::linalg::solve(xt::transpose(v), xt::transpose(gvt)));
664
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 4 return xt::stack(xt::xtuple(dr, ds, dt), 2);
665 2 }
666
667 } // namespace oiseau::dg::nodal::utils
668