GCC Code Coverage Report

Directory: src/oiseau/
File: src/oiseau/dg/nodal/ref_tetrahedron.cpp
Date: 2025-05-24 01:28:39
Exec Total Coverage
Lines: 232 232 100.0%
Functions: 13 13 100.0%
Branches: 293 530 55.3%
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_tetrahedron.hpp"
8
9 #include <cmath>
10 #include <cstddef>
11 #include <numbers>
12 #include <xtensor-blas/xlinalg.hpp>
13 #include <xtensor/containers/xtensor.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 RefTetrahedron::RefTetrahedron(unsigned order) : RefElement(order) {
30 16 this->m_np = ((order + 1) * (order + 2) * (order + 3)) / 6;
31 16 this->m_nfp = ((order + 1) * (order + 2)) / 2;
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_xyz_to_rst(detail::generate_tetrahedron_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 621 xt::xarray<double> RefTetrahedron::basis_function(const xt::xarray<double> &abc, int i, int j,
39 int k) const {
40
2/4
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 xt::xarray<double> a = xt::col(abc, 0);
41
2/4
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 xt::xarray<double> b = xt::col(abc, 1);
42
2/4
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 xt::xarray<double> c = xt::col(abc, 2);
43
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
44
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 auto h2 = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
45
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 auto h3 = oiseau::utils::jacobi_p(k, 2.0 * (i + j) + 2.0, 0.0, c);
46
1/2
✓ Branch 10 taken 621 times.
✗ Branch 11 not taken.
 1242 return 2 * std::numbers::sqrt2 * h1 * h2 * xt::pow(1 - b, i) * h3 * xt::pow(1 - c, i + j);
47 621 }
48
49 621 xt::xarray<double> RefTetrahedron::grad_basis_function(const xt::xarray<double> &abc, int i, int j,
50 int k) const {
51
2/4
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 xt::xarray<double> a = xt::col(abc, 0);
52
2/4
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 xt::xarray<double> b = xt::col(abc, 1);
53
2/4
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 xt::xarray<double> c = xt::col(abc, 2);
54
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 xt::xarray<double> fa = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
55
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 xt::xarray<double> dfa = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, a);
56
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 xt::xarray<double> gb = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
57
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 xt::xarray<double> dgb = oiseau::utils::grad_jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
58
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 xt::xarray<double> hc = oiseau::utils::jacobi_p(k, 2.0 * (i + j) + 2.0, 0.0, c);
59
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 xt::xarray<double> dhc = oiseau::utils::grad_jacobi_p(k, 2.0 * (i + j) + 2.0, 0.0, c);
60
61
4/8
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 621 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 621 times.
✗ Branch 11 not taken.
 621 xt::xarray<double> v3dr, v3ds, v3dt, tmp;
62
63
1/2
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
 621 v3dr = dfa * (gb * hc);
64
3/4
✓ Branch 0 taken 380 times.
✓ Branch 1 taken 241 times.
✓ Branch 7 taken 380 times.
✗ Branch 8 not taken.
 621 if (i > 0) v3dr = v3dr * xt::pow(0.5 * (1 - b), i - 1);
65
3/4
✓ Branch 0 taken 545 times.
✓ Branch 1 taken 76 times.
✓ Branch 7 taken 545 times.
✗ Branch 8 not taken.
 621 if (i + j > 0) v3dr = v3dr * xt::pow(0.5 * (1 - c), i + j - 1);
66
67
1/2
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
 621 v3ds = 0.5 * (1 + a) * v3dr;
68
1/2
✓ Branch 5 taken 621 times.
✗ Branch 6 not taken.
 621 tmp = dgb * xt::pow(0.5 * (1 - b), i);
69
3/4
✓ Branch 0 taken 380 times.
✓ Branch 1 taken 241 times.
✓ Branch 9 taken 380 times.
✗ Branch 10 not taken.
 621 if (i > 0) tmp = tmp + (-0.5 * i) * gb * xt::pow(0.5 * (1 - b), i - 1);
70
3/4
✓ Branch 0 taken 545 times.
✓ Branch 1 taken 76 times.
✓ Branch 7 taken 545 times.
✗ Branch 8 not taken.
 621 if (i + j > 0) tmp = tmp * xt::pow(0.5 * (1 - c), i + j - 1);
71
1/2
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
 621 tmp = fa * (tmp * hc);
72
1/2
✓ Branch 2 taken 621 times.
✗ Branch 3 not taken.
 621 v3ds = v3ds + tmp;
73
74
1/2
✓ Branch 8 taken 621 times.
✗ Branch 9 not taken.
 621 v3dt = 0.5 * (1 + a) * v3dr + 0.5 * (1 + b) * tmp;
75
1/2
✓ Branch 5 taken 621 times.
✗ Branch 6 not taken.
 621 tmp = dhc * xt::pow(0.5 * (1 - c), i + j);
76
3/4
✓ Branch 0 taken 545 times.
✓ Branch 1 taken 76 times.
✓ Branch 9 taken 545 times.
✗ Branch 10 not taken.
 621 if (i + j > 0) tmp = tmp - 0.5 * (i + j) * hc * xt::pow(0.5 * (1 - c), (i + j - 1));
77
1/2
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
 621 tmp = fa * gb * tmp;
78
1/2
✓ Branch 5 taken 621 times.
✗ Branch 6 not taken.
 621 tmp = tmp * xt::pow(0.5 * (1 - b), i);
79
1/2
✓ Branch 2 taken 621 times.
✗ Branch 3 not taken.
 621 v3dt = v3dt + tmp;
80
81
1/2
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
 621 v3dr = v3dr * std::pow(2, 2 * i + j + 1.5);
82
1/2
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
 621 v3ds = v3ds * std::pow(2, 2 * i + j + 1.5);
83
1/2
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
 621 v3dt = v3dt * std::pow(2, 2 * i + j + 1.5);
84
85
2/4
✓ Branch 2 taken 621 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 621 times.
✗ Branch 6 not taken.
 1242 return xt::stack(xt::xtuple(v3dr, v3ds, v3dt), 1);
86 621 }
87
88 16 xt::xarray<double> RefTetrahedron::vandermonde(const xt::xarray<double> &rst) const {
89
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto abc = detail::rst_to_abc(rst);
90 16 const std::size_t n_points = rst.shape()[0];
91 16 const std::size_t n_basis = this->m_np; // (order+1)(order+2)*(order+3)/6
92 //
93
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 xt::xarray<double> output(xt::dynamic_shape<std::size_t>({n_points, n_basis}));
94
95 16 std::size_t index = 0;
96
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (unsigned i = 0; i <= this->m_order; ++i) {
97
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 76 times.
 317 for (unsigned j = 0; j <= this->m_order - i; ++j) {
98
2/2
✓ Branch 0 taken 621 times.
✓ Branch 1 taken 241 times.
 862 for (unsigned k = 0; k <= this->m_order - i - j; ++k, ++index) {
99
3/6
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 621 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 621 times.
✗ Branch 8 not taken.
 621 xt::col(output, index) = this->basis_function(abc, i, j, k);
100 }
101 }
102 }
103 16 return output;
104 16 }
105
106 16 xt::xarray<double> RefTetrahedron::grad_vandermonde(const xt::xarray<double> &rst) const {
107
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto abc = detail::rst_to_abc(rst);
108 16 const std::size_t n_points = rst.shape()[0];
109 16 const std::size_t n_basis = this->m_np; // (order+1)(order+2)*(order+3)/6
110 //
111
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 xt::xarray<double> output(xt::dynamic_shape<std::size_t>({n_points, n_basis, 3}));
112
113 16 std::size_t index = 0;
114
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (unsigned i = 0; i <= this->m_order; ++i) {
115
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 76 times.
 317 for (unsigned j = 0; j <= this->m_order - i; ++j) {
116
2/2
✓ Branch 0 taken 621 times.
✓ Branch 1 taken 241 times.
 862 for (unsigned k = 0; k <= this->m_order - i - j; ++k, ++index) {
117
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 auto tmp = this->grad_basis_function(abc, i, j, k);
118
2/4
✓ Branch 3 taken 621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 621 times.
✗ Branch 7 not taken.
 621 xt::view(output, xt::all(), index, xt::all()) = tmp;
119 621 }
120 }
121 }
122 16 return output;
123 16 }
124 //
125 16 xt::xarray<double> RefTetrahedron::grad_operator(const xt::xarray<double> &v,
126 const xt::xarray<double> &gv) const {
127
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
128
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvs = xt::view(gv, xt::all(), xt::all(), 1);
129
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvt = xt::view(gv, xt::all(), xt::all(), 2);
130
131 16 const auto vt = xt::transpose(v);
132 16 const auto dvt_dr = xt::transpose(gvr);
133 16 const auto dvt_ds = xt::transpose(gvs);
134 16 const auto dvt_dt = xt::transpose(gvt);
135
136
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto dr = xt::transpose(xt::linalg::solve(vt, dvt_dr));
137
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto ds = xt::transpose(xt::linalg::solve(vt, dvt_ds));
138
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto dt = xt::transpose(xt::linalg::solve(vt, dvt_dt));
139
140
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, dt), 2);
141 16 }
142
143 } // namespace oiseau::dg::nodal
144
145 namespace oiseau::dg::nodal::detail {
146
147 32 xt::xarray<double> rst_to_abc(const xt::xarray<double> &rst) {
148
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto abc = xt::zeros_like(rst);
149
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto &&r = xt::col(rst, 0);
150
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto &&s = xt::col(rst, 1);
151
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto &&t = xt::col(rst, 2);
152
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto &&a = xt::col(abc, 0);
153
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto &&b = xt::col(abc, 1);
154
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto &&c = xt::col(abc, 2);
155
2/2
✓ Branch 2 taken 1242 times.
✓ Branch 3 taken 32 times.
 1274 for (unsigned n = 0; n < rst.shape()[0]; ++n) {
156
4/6
✓ Branch 1 taken 1242 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1242 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1136 times.
✓ Branch 7 taken 106 times.
 1242 if ((s[n] + t[n]) != 0) {
157
4/8
✓ Branch 1 taken 1136 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1136 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1136 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1136 times.
✗ Branch 11 not taken.
 1136 a[n] = 2 * ((1 + r[n]) / (-s[n] - t[n])) - 1;
158 } else {
159
1/2
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
 106 a[n] = -1;
160 }
161
3/4
✓ Branch 1 taken 1242 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1210 times.
✓ Branch 4 taken 32 times.
 1242 if (t[n] != 1) {
162
3/6
✓ Branch 1 taken 1210 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1210 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1210 times.
✗ Branch 8 not taken.
 1210 b[n] = 2 * (1 + s[n]) / (1 - t[n]) - 1;
163 } else {
164
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 b[n] = -1;
165 }
166 }
167
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 c = t;
168 32 return abc;
169 32 }
170
171 192 xt::xarray<double> eval_warp(int p, const xt::xarray<double> &xnodes,
172 const xt::xarray<double> &xout) {
173
1/2
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
 192 xt::xarray<double> warp = xt::zeros_like(xout);
174
1/2
✓ Branch 2 taken 192 times.
✗ Branch 3 not taken.
 192 xt::xarray<double> xeq = xt::zeros<double>({p + 1});
175
3/4
✓ Branch 1 taken 912 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 912 times.
✓ Branch 4 taken 192 times.
 1104 for (int i = 0; i <= p; i++) xeq(i) = -1.0 + 2.0 * (p - i) / p;
176
2/2
✓ Branch 0 taken 912 times.
✓ Branch 1 taken 192 times.
 1104 for (int i = 0; i <= p; i++) {
177
4/8
✓ Branch 1 taken 912 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 912 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 912 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 912 times.
✗ Branch 12 not taken.
 912 xt::xarray<double> d = xt::ones_like(xout) * (xnodes(i) - xeq(i));
178
2/2
✓ Branch 0 taken 3048 times.
✓ Branch 1 taken 912 times.
 3960 for (int j = 1; j < p; j++) {
179
6/10
✓ Branch 0 taken 2520 times.
✓ Branch 1 taken 528 times.
✓ Branch 3 taken 2520 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2520 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2520 times.
✗ Branch 10 not taken.
✓ Branch 15 taken 2520 times.
✗ Branch 16 not taken.
 3048 if (i != j) d = d * (xout - xeq(j)) / (xeq(i) - xeq(j));
180 }
181
5/8
✓ Branch 0 taken 720 times.
✓ Branch 1 taken 192 times.
✓ Branch 3 taken 720 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 720 times.
✗ Branch 7 not taken.
✓ Branch 11 taken 720 times.
✗ Branch 12 not taken.
 912 if (i != 0) d = -d / (xeq(i) - xeq(0));
182
5/8
✓ Branch 0 taken 720 times.
✓ Branch 1 taken 192 times.
✓ Branch 3 taken 720 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 720 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 720 times.
✗ Branch 11 not taken.
 912 if (i != p) d = d / (xeq(i) - xeq(p));
183
1/2
✓ Branch 1 taken 912 times.
✗ Branch 2 not taken.
 912 warp += d;
184 912 }
185 192 return warp;
186 192 }
187
188 64 xt::xarray<double> eval_shift(int p, double pval, const xt::xarray<double> &l1,
189 const xt::xarray<double> &l2, const xt::xarray<double> &l3) {
190
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 auto gaussX = -oiseau::dg::nodal::utils::jacobi_gl(p, 0.0, 0.0);
191 64 auto blend1 = l2 * l3;
192 64 auto blend2 = l1 * l3;
193 64 auto blend3 = l1 * l2;
194
3/6
✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 64 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
 64 auto warpfactor1 = 4 * eval_warp(p, gaussX, l3 - 1 * l2);
195
3/6
✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 64 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
 64 auto warpfactor2 = 4 * eval_warp(p, gaussX, l1 - 1 * l3);
196
3/6
✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 64 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
 64 auto warpfactor3 = 4 * eval_warp(p, gaussX, l2 - 1 * l1);
197 64 auto warp1 = blend1 * warpfactor1 * (1 + xt::pow(pval * l1, 2));
198 64 auto warp2 = blend2 * warpfactor2 * (1 + xt::pow(pval * l2, 2));
199 64 auto warp3 = blend3 * warpfactor3 * (1 + xt::pow(pval * l3, 2));
200 64 auto pi = xt::numeric_constants<double>::PI;
201 64 auto dx = 1 * warp1 + cos(2.0 * pi / 3.0) * warp2 + cos(4.0 * pi / 3.0) * warp3;
202 64 auto dy = 0 * warp1 + sin(2.0 * pi / 3.0) * warp2 + sin(4.0 * pi / 3.0) * warp3;
203
2/4
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
 128 return xt::stack(xt::xtuple(dx, dy), 1);
204 64 }
205
206 64 xt::xarray<double> warp_shift_face_3d(int p, double pval, const xt::xarray<double> &l2,
207 const xt::xarray<double> &l3, const xt::xarray<double> &l4) {
208 64 return eval_shift(p, pval, l2, l3, l4);
209 }
210
211 16 xt::xarray<double> generate_tetrahedron_equidistant_nodes(unsigned n) {
212 16 auto n_p = ((n + 1) * (n + 2) * (n + 3)) / 6;
213
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 auto shape = xt::xarray<double>::shape_type{n_p, 3};
214
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 xt::xarray<double> output = xt::zeros<double>(shape);
215 16 std::size_t sk = 0;
216
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (std::size_t i = 0; i < n + 1; ++i) {
217
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 76 times.
 317 for (std::size_t j = 0; j < n + 1 - i; ++j) {
218
2/2
✓ Branch 0 taken 621 times.
✓ Branch 1 taken 241 times.
 862 for (std::size_t k = 0; k < n + 1 - i - j; ++k, ++sk) {
219
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 output(sk, 0) = -1 + k * 2.0 / n;
220
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 output(sk, 1) = -1 + j * 2.0 / n;
221
1/2
✓ Branch 1 taken 621 times.
✗ Branch 2 not taken.
 621 output(sk, 2) = -1 + i * 2.0 / n;
222 }
223 }
224 }
225 16 return output;
226 16 }
227
228 16 xt::xarray<double> generate_tetrahedron_nodes(unsigned p) {
229 xt::xarray<double> alphastore = {
230 0, 0, 0, 0.1002, 1.1332, 1.5608, 1.3413, 1.2577,
231 1.1603, 1.10153, 0.6080, 0.4523, 0.8856, 0.8717, 0.9655,
232
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 };
233
2/4
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 double alpha = p > 14 ? 1 : alphastore[p];
234 16 double tol = 1e-10;
235
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto rst = generate_tetrahedron_equidistant_nodes(p);
236
3/6
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 auto r = xt::col(rst, 0), s = xt::col(rst, 1), t = xt::col(rst, 2);
237
3/6
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 10 taken 16 times.
✗ Branch 11 not taken.
✓ Branch 20 taken 16 times.
✗ Branch 21 not taken.
 16 xt::xarray<double> l1 = (1 + t) / 2, l2 = (1 + s) / 2, l3 = -(1 + r + s + t) / 2,
238
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 l4 = (1 + r) / 2;
239
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v1 = {-1.0, -1.0 / std::numbers::sqrt3, -1.0 / std::sqrt(6)};
240
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v2 = {1.0, -1.0 / std::numbers::sqrt3, -1.0 / std::sqrt(6)};
241
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v3 = {0.0, 2.0 / std::numbers::sqrt3, -1.0 / std::sqrt(6)};
242
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v4 = {0.0, 0.0, 3.0 / std::sqrt(6)};
243
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xtensor<double, 2> t1 = xt::zeros<double>({4, 3});
244
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xtensor<double, 2> t2 = xt::zeros<double>({4, 3});
245
2/4
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
 16 xt::row(t1, 0) = v2 - 1 * v1;
246
2/4
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
 16 xt::row(t1, 1) = v2 - 1 * v1;
247
2/4
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
 16 xt::row(t1, 2) = v3 - 1 * v2;
248
2/4
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
 16 xt::row(t1, 3) = v3 - 1 * v1;
249
2/4
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 xt::row(t2, 0) = v3 - 0.5 * (v1 + v2);
250
2/4
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 xt::row(t2, 1) = v4 - 0.5 * (v1 + v2);
251
2/4
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 xt::row(t2, 2) = v4 - 0.5 * (v2 + v3);
252
2/4
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 xt::row(t2, 3) = v4 - 0.5 * (v1 + v3);
253
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 16 times.
 80 for (int i = 0; i < 4; ++i) {
254
4/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 10 taken 64 times.
✗ Branch 11 not taken.
 64 xt::row(t1, i) /= xt::linalg::norm(xt::row(t1, i));
255
4/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 10 taken 64 times.
✗ Branch 11 not taken.
 64 xt::row(t2, i) /= xt::linalg::norm(xt::row(t2, i));
256 }
257
3/6
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 32 auto xyz = xt::linalg::outer(l3, v1) + xt::linalg::outer(l4, v2) + xt::linalg::outer(l2, v3) +
258
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 48 xt::linalg::outer(l1, v4);
259
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 xt::xarray<double> shift = xt::zeros<double>(xyz.shape());
260
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 16 times.
 80 for (int face = 0; face < 4; ++face) {
261
4/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 10 taken 64 times.
✗ Branch 11 not taken.
 64 xt::xarray<double> la, lb, lc, ld;
262
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 48 times.
 64 if (face == 0) {
263
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 la = l1;
264
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lb = l2;
265
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lc = l3;
266
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ld = l4;
267
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 32 times.
 48 } else if (face == 1) {
268
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 la = l2;
269
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lb = l1;
270
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lc = l3;
271
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ld = l4;
272
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 16 times.
 32 } else if (face == 2) {
273
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 la = l3;
274
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lb = l1;
275
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lc = l4;
276
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ld = l2;
277
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 } else if (face == 3) {
278
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 la = l4;
279
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lb = l1;
280
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 lc = l3;
281
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ld = l2;
282 }
283
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 auto warp = warp_shift_face_3d(p, alpha, lb, lc, ld);
284
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 auto warp1 = xt::col(warp, 0);
285
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 auto warp2 = xt::col(warp, 1);
286
1/2
✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
 64 xt::xarray<double> blend = lb * lc * ld;
287
1/2
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
 64 xt::xarray<double> denom = (lb + 0.5 * la) * (lc + 0.5 * la) * (ld + 0.5 * la);
288
2/4
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
 64 auto ids = xt::flatten_indices(xt::argwhere(denom > tol));
289
4/8
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 64 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 64 times.
✗ Branch 14 not taken.
 256 xt::view(blend, xt::keep(ids)) = (1 + xt::pow(alpha * xt::view(la, xt::keep(ids)), 2)) *
290
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 256 xt::view(blend, xt::keep(ids)) /
291
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 256 xt::view(denom, xt::keep(ids));
292
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
 128 shift = shift + xt::linalg::outer(blend * warp1, xt::row(t1, face)) +
293
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
 192 xt::linalg::outer(blend * warp2, xt::row(t2, face));
294 64 auto mask = (la < tol) && (((lb > tol) + (lc > tol) + (ld > tol)) < 3);
295
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 64 ids = xt::flatten_indices(xt::argwhere(mask));
296
2/4
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
 128 xt::view(shift, xt::keep(ids), xt::all()) =
297
4/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 10 taken 64 times.
✗ Branch 11 not taken.
 128 xt::linalg::outer(xt::view(warp1, xt::keep(ids)), xt::row(t1, face)) +
298
5/10
✓ 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 10 taken 64 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 64 times.
✗ Branch 15 not taken.
 256 xt::linalg::outer(xt::view(warp2, xt::keep(ids)), xt::row(t2, face));
299 64 }
300 16 xyz = xyz + shift;
301
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 return xyz;
302 16 }
303
304 16 xt::xarray<double> equilateral_xyz_to_rst(const xt::xarray<double> &coords) {
305 16 const double sqrt3 = std::numbers::sqrt3;
306 16 const double sqrt6 = std::sqrt(6.0);
307
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v1 = {-1, -1 / sqrt3, -1 / sqrt6};
308
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v2 = {1, -1 / sqrt3, -1 / sqrt6};
309
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v3 = {0, 2 / sqrt3, -1 / sqrt6};
310
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xtensor<double, 1> v4 = {0, 0 / sqrt3, 3 / sqrt6};
311 16 auto rhs = xt::transpose(coords - 0.5 * (v2 + v3 + v4 - v1));
312
1/2
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 auto lhs = 0.5 * xt::stack(xt::xtuple(v2 - v1, v3 - v1, v4 - v1), 1);
313
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto rst = xt::linalg::solve(lhs, rhs);
314
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 32 return xt::transpose(rst);
315 16 }
316
317 } // namespace oiseau::dg::nodal::detail
318