GCC Code Coverage Report

Directory: src/oiseau/
File: src/oiseau/dg/nodal/ref_quadrilateral.cpp
Date: 2025-05-24 01:28:39
Exec Total Coverage
Lines: 61 63 96.8%
Functions: 7 7 100.0%
Branches: 50 92 54.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_quadrilateral.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 16 RefQuadrilateral::RefQuadrilateral(unsigned order) : RefElement(order) {
25 16 this->m_np = (order + 1) * (order + 1);
26 16 this->m_nfp = order + 1;
27
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_r = detail::generate_quadrilateral_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 = grad_operator(this->m_v, this->m_gv);
31 16 }
32
33 406 xt::xarray<double> RefQuadrilateral::basis_function(const xt::xarray<double> &rs, int i, int j) {
34
2/4
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 406 times.
✗ Branch 5 not taken.
 406 xt::xarray<double> r = xt::col(rs, 0);
35
2/4
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 406 times.
✗ Branch 5 not taken.
 406 xt::xarray<double> s = xt::col(rs, 1);
36
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, r);
37
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto h2 = oiseau::utils::jacobi_p(j, 0.0, 0.0, s);
38
1/2
✓ Branch 2 taken 406 times.
✗ Branch 3 not taken.
 1218 return h1 * h2;
39 406 }
40 //
41 406 xt::xarray<double> RefQuadrilateral::grad_basis_function(const xt::xarray<double> &rs, int i,
42 int j) {
43
2/4
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 406 times.
✗ Branch 5 not taken.
 406 xt::xarray<double> r = xt::col(rs, 0);
44
2/4
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 406 times.
✗ Branch 5 not taken.
 406 xt::xarray<double> s = xt::col(rs, 1);
45
46
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto pi_r = oiseau::utils::jacobi_p(i, 0.0, 0.0, r); // P_i(r)
47
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto dpi_r = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, r); // dP_i/dr
48
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto pj_s = oiseau::utils::jacobi_p(j, 0.0, 0.0, s); // P_j(s)
49
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto dpj_s = oiseau::utils::grad_jacobi_p(j, 0.0, 0.0, s); // dP_j/ds
50
51 406 auto dphidr = dpi_r * pj_s; // ∂φ/∂r
52 406 auto dphids = pi_r * dpj_s; // ∂φ/∂s
53
54
2/4
✓ Branch 2 taken 406 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 406 times.
✗ Branch 6 not taken.
 812 return xt::stack(xt::xtuple(dphidr, dphids), 1); // shape: [N, 2]
55 406 }
56
57 16 xt::xarray<double> RefQuadrilateral::vandermonde(const xt::xarray<double> &rs) const {
58 16 const std::size_t n_points = rs.shape()[0];
59 16 const std::size_t n_basis = this->m_np; // (order+1)*(order+1)
60
61
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}));
62 16 int index = 0;
63
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (unsigned i = 0; i <= this->m_order; ++i) {
64
2/2
✓ Branch 0 taken 406 times.
✓ Branch 1 taken 76 times.
 482 for (unsigned j = 0; j <= this->m_order; ++j, ++index) {
65
3/6
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 406 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 406 times.
✗ Branch 8 not taken.
 406 xt::col(output, index) = this->basis_function(rs, i, j);
66 }
67 }
68 16 return output;
69 ✗ }
70
71 16 xt::xarray<double> RefQuadrilateral::grad_vandermonde(const xt::xarray<double> &rs) const {
72 16 const std::size_t n_points = rs.shape()[0];
73 16 const std::size_t n_basis = this->m_np; // (order+1)*(order+1)
74
75
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, 2}));
76
77 16 int index = 0;
78
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (unsigned i = 0; i <= this->m_order; ++i) {
79
2/2
✓ Branch 0 taken 406 times.
✓ Branch 1 taken 76 times.
 482 for (unsigned j = 0; j <= this->m_order; ++j, ++index) {
80
1/2
✓ Branch 1 taken 406 times.
✗ Branch 2 not taken.
 406 auto tmp = this->grad_basis_function(rs, i, j);
81
2/4
✓ Branch 3 taken 406 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 406 times.
✗ Branch 7 not taken.
 406 xt::view(output, xt::all(), index, xt::all()) = tmp;
82 406 }
83 }
84 16 return output;
85 ✗ }
86
87 16 xt::xarray<double> RefQuadrilateral::grad_operator(const xt::xarray<double> &v,
88 const xt::xarray<double> &gv) const {
89
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
90
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvs = xt::view(gv, xt::all(), xt::all(), 1);
91
92 16 const auto vt = xt::transpose(v);
93 16 const auto dvt_dr = xt::transpose(gvr);
94 16 const auto dvt_ds = xt::transpose(gvs);
95
96
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto dr = xt::transpose(xt::linalg::solve(vt, dvt_dr));
97
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto ds = xt::transpose(xt::linalg::solve(vt, dvt_ds));
98
99
2/4
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 32 return xt::stack(xt::xtuple(dr, ds), 2);
100 16 }
101
102 } // namespace oiseau::dg::nodal
103
104 namespace oiseau::dg::nodal::detail {
105
106 16 xt::xarray<double> generate_quadrilateral_nodes(unsigned order) {
107
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 xt::xarray<double> r1d = utils::jacobi_gl(order, 0.0, 0.0);
108 16 auto [R, S] = xt::meshgrid(r1d, r1d);
109
4/8
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 16 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 16 times.
✗ Branch 12 not taken.
 32 return xt::stack(xt::xtuple(xt::flatten(R), xt::flatten(S)), 1);
110 16 }
111 //
112
113 } // namespace oiseau::dg::nodal::detail
114