GCC Code Coverage Report

Directory: src/oiseau/
File: src/oiseau/dg/nodal/ref_line.cpp
Date: 2025-05-24 01:28:39
Exec Total Coverage
Lines: 39 41 95.1%
Functions: 7 7 100.0%
Branches: 27 48 56.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_line.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/utils/integration.hpp"
20 #include "oiseau/utils/math.hpp"
21
22 namespace oiseau::dg::nodal {
23
24 17 RefLine::RefLine(unsigned order) : RefElement(order) {
25 16 this->m_np = order + 1;
26 16 this->m_nfp = 1;
27
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_r = detail::generate_line_nodes(this->m_order);
28
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_v = this->vandermonde(this->m_r);
29
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_gv = this->grad_vandermonde(this->m_r);
30
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_d = this->grad_operator(this->m_v, this->m_gv);
31 16 }
32
33 76 xt::xarray<double> RefLine::basis_function(const xt::xarray<double>& r, int i) {
34 76 return oiseau::utils::jacobi_p(i, 0.0, 0.0, r);
35 }
36
37 76 xt::xarray<double> RefLine::grad_basis_function(const xt::xarray<double>& r, int i) {
38 76 return oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, r);
39 }
40
41 16 xt::xarray<double> RefLine::vandermonde(const xt::xarray<double>& r) const {
42 16 const auto n_points = r.shape()[0];
43 16 const auto n_basis = this->m_np; // order + 1
44
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}));
45
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (std::size_t i = 0; i <= this->m_order; ++i) {
46
3/6
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 76 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 76 times.
✗ Branch 9 not taken.
 76 xt::view(output, xt::all(), i) = this->basis_function(r, i);
47 }
48 16 return output;
49 ✗ }
50
51 16 xt::xarray<double> RefLine::grad_vandermonde(const xt::xarray<double>& r) const {
52 16 const auto n_points = r.shape()[0];
53 16 const auto n_basis = this->m_np; // order + 1
54
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}));
55
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (size_t i = 0; i <= this->m_order; ++i) {
56
3/6
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 76 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 76 times.
✗ Branch 9 not taken.
 76 xt::view(output, xt::all(), i) = grad_basis_function(r, i);
57 }
58 16 return output;
59 ✗ }
60
61 16 xt::xarray<double> RefLine::grad_operator(const xt::xarray<double>& v,
62 const xt::xarray<double>& gv) const {
63
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
64
65 16 const auto vt = xt::transpose(v);
66 16 const auto dvt_dr = xt::transpose(gvr);
67
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto dr = xt::transpose(xt::linalg::solve(vt, dvt_dr));
68
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 return dr;
69 16 }
70
71 } // namespace oiseau::dg::nodal
72
73 namespace oiseau::dg::nodal::detail {
74
75 16 xt::xarray<double> generate_line_nodes(unsigned order) {
76
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 14 times.
 16 if (order == 1) {
77
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return xt::xarray<double>{-1.0, 1.0};
78 }
79
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 auto [x, _] = oiseau::utils::integration::jacobi_gq(order - 2, 1.0, 1.0);
80
2/4
✓ Branch 5 taken 14 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 14 times.
✗ Branch 9 not taken.
 14 return xt::concatenate(xt::xtuple(-xt::ones<double>({1}), x, xt::ones<double>({1})));
81 14 }
82 } // namespace oiseau::dg::nodal::detail
83