GCC Code Coverage Report

Directory: src/oiseau/
File: src/oiseau/dg/nodal/ref_hexahedron.cpp
Date: 2025-05-24 01:28:39
Exec Total Coverage
Lines: 72 74 97.3%
Functions: 7 7 100.0%
Branches: 64 116 55.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_hexahedron.hpp"
8
9 #include <cstddef>
10 #include <xtensor-blas/xlinalg.hpp>
11 #include <xtensor/core/xshape.hpp>
12 #include <xtensor/core/xtensor_forward.hpp>
13 #include <xtensor/generators/xbuilder.hpp>
14 #include <xtensor/misc/xmanipulation.hpp>
15 #include <xtensor/views/xslice.hpp>
16 #include <xtensor/views/xview.hpp>
17
18 #include "oiseau/dg/nodal/ref_element.hpp"
19 #include "oiseau/dg/nodal/utils.hpp"
20 #include "oiseau/utils/math.hpp"
21
22 namespace oiseau::dg::nodal {
23
24 14 RefHexahedron::RefHexahedron(unsigned order) : RefElement(order) {
25 14 this->m_np = (order + 1) * (order + 1) * (order + 1);
26 14 this->m_nfp = (order + 1) * (order + 1);
27
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 this->m_r = detail::generate_hexahedron_nodes(this->m_order);
28
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 this->m_v = this->vandermonde(this->m_r);
29
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 this->m_gv = this->grad_vandermonde(this->m_r);
30
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 this->m_d = this->grad_operator(this->m_v, this->m_gv);
31 14 }
32
33 1312 xt::xarray<double> RefHexahedron::basis_function(const xt::xarray<double> &rst, int i, int j,
34 int k) const {
35
2/4
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
 1312 xt::xarray<double> r = xt::col(rst, 0);
36
2/4
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
 1312 xt::xarray<double> s = xt::col(rst, 1);
37
2/4
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
 1312 xt::xarray<double> t = xt::col(rst, 2);
38
39
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, r); // P_i(r)
40
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto h2 = oiseau::utils::jacobi_p(j, 0.0, 0.0, s); // P_j(s)
41
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto h3 = oiseau::utils::jacobi_p(k, 0.0, 0.0, t); // P_k(t)
42
43
1/2
✓ Branch 3 taken 1312 times.
✗ Branch 4 not taken.
 2624 return h1 * h2 * h3;
44 1312 }
45
46 1312 xt::xarray<double> RefHexahedron::grad_basis_function(const xt::xarray<double> &rst, int i, int j,
47 int k) const {
48
2/4
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
 1312 xt::xarray<double> r = xt::col(rst, 0);
49
2/4
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
 1312 xt::xarray<double> s = xt::col(rst, 1);
50
2/4
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
 1312 xt::xarray<double> t = xt::col(rst, 2);
51
52
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto Pi_r = oiseau::utils::jacobi_p(i, 0.0, 0.0, r); // P_i(r)
53
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto Pj_s = oiseau::utils::jacobi_p(j, 0.0, 0.0, s); // P_j(s)
54
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto Pk_t = oiseau::utils::jacobi_p(k, 0.0, 0.0, t); // P_k(t)
55
56
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto dPi_r = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, r); // dP_i/dr
57
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto dPj_s = oiseau::utils::grad_jacobi_p(j, 0.0, 0.0, s); // dP_j/ds
58
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto dPk_t = oiseau::utils::grad_jacobi_p(k, 0.0, 0.0, t); // dP_k/dt
59
60 1312 auto dphidr = dPi_r * Pj_s * Pk_t;
61 1312 auto dphids = Pi_r * dPj_s * Pk_t;
62 1312 auto dphidt = Pi_r * Pj_s * dPk_t;
63
64
2/4
✓ Branch 2 taken 1312 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1312 times.
✗ Branch 6 not taken.
 2624 return xt::stack(xt::xtuple(dphidr, dphids, dphidt), 1); // shape: [N, 3]
65 1312 }
66
67 14 xt::xarray<double> RefHexahedron::vandermonde(const xt::xarray<double> &rst) const {
68 14 const std::size_t n_points = rst.shape()[0];
69 14 const std::size_t n_basis = this->m_np; // (order+1)*(order+1)*(order+1)
70
71
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 14 xt::xarray<double> output(xt::dynamic_shape<std::size_t>({n_points, n_basis}));
72 14 int index = 0;
73
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 14 times.
 72 for (unsigned i = 0; i <= this->m_order; ++i) {
74
2/2
✓ Branch 0 taken 266 times.
✓ Branch 1 taken 58 times.
 324 for (unsigned j = 0; j <= this->m_order; ++j) {
75
2/2
✓ Branch 0 taken 1312 times.
✓ Branch 1 taken 266 times.
 1578 for (unsigned k = 0; k <= this->m_order; ++k, ++index) {
76
3/6
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1312 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1312 times.
✗ Branch 8 not taken.
 1312 xt::col(output, index) = this->basis_function(rst, i, j, k);
77 }
78 }
79 }
80 14 return output;
81 ✗ }
82
83 14 xt::xarray<double> RefHexahedron::grad_vandermonde(const xt::xarray<double> &rst) const {
84 14 const std::size_t n_points = rst.shape()[0];
85 14 const std::size_t n_basis = this->m_np; // (order+1)*(order+1)*(order+1)
86
87
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 14 xt::xarray<double> output(xt::dynamic_shape<std::size_t>({n_points, n_basis, 3}));
88 14 std::size_t index = 0;
89
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 14 times.
 72 for (unsigned i = 0; i <= this->m_order; ++i) {
90
2/2
✓ Branch 0 taken 266 times.
✓ Branch 1 taken 58 times.
 324 for (unsigned j = 0; j <= this->m_order; ++j) {
91
2/2
✓ Branch 0 taken 1312 times.
✓ Branch 1 taken 266 times.
 1578 for (unsigned k = 0; k <= this->m_order; ++k, ++index) {
92
1/2
✓ Branch 1 taken 1312 times.
✗ Branch 2 not taken.
 1312 auto tmp = this->grad_basis_function(rst, i, j, k);
93
2/4
✓ Branch 3 taken 1312 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1312 times.
✗ Branch 7 not taken.
 1312 xt::view(output, xt::all(), index, xt::all()) = tmp;
94 1312 }
95 }
96 }
97 14 return output;
98 ✗ }
99
100 14 xt::xarray<double> RefHexahedron::grad_operator(const xt::xarray<double> &v,
101 const xt::xarray<double> &gv) const {
102
1/2
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
 14 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
103
1/2
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
 14 auto gvs = xt::view(gv, xt::all(), xt::all(), 1);
104
1/2
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
 14 auto gvt = xt::view(gv, xt::all(), xt::all(), 2);
105
106 14 const auto vt = xt::transpose(v);
107 14 const auto dvt_dr = xt::transpose(gvr);
108 14 const auto dvt_ds = xt::transpose(gvs);
109 14 const auto dvt_dt = xt::transpose(gvt);
110
111
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 const auto dr = xt::transpose(xt::linalg::solve(vt, dvt_dr));
112
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 const auto ds = xt::transpose(xt::linalg::solve(vt, dvt_ds));
113
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 const auto dt = xt::transpose(xt::linalg::solve(vt, dvt_dt));
114
115
2/4
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 14 times.
✗ Branch 6 not taken.
 28 return xt::stack(xt::xtuple(dr, ds, dt), 2);
116 14 }
117
118 } // namespace oiseau::dg::nodal
119
120 namespace oiseau::dg::nodal::detail {
121 14 xt::xarray<double> generate_hexahedron_nodes(unsigned order) {
122
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 xt::xarray<double> r1d = utils::jacobi_gl(order, 0.0, 0.0);
123 14 auto [R, S, T] = xt::meshgrid(r1d, r1d, r1d);
124
5/10
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 14 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 14 times.
✗ Branch 15 not taken.
 28 return xt::stack(xt::xtuple(xt::flatten(R), xt::flatten(S), xt::flatten(T)), 1);
125 14 }
126
127 } // namespace oiseau::dg::nodal::detail
128