GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/LinearOperator/DiscreteLocalOperator.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 42 44 95.5%
Functions: 4 4 100.0%
Branches: 32 52 61.5%
Line Branch Exec Source
1 // This file is part of Bembel, the higher order C++ boundary element library.
2 //
3 // Copyright (C) 2022 see <http://www.bembel.eu>
4 //
5 // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
6 // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
7 // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
8 // source code is subject to the GNU General Public License version 3 and
9 // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
10 // information.
11
12 #ifndef BEMBEL_SRC_LINEAROPERATOR_DISCRETELOCALOPERATOR_HPP_
13 #define BEMBEL_SRC_LINEAROPERATOR_DISCRETELOCALOPERATOR_HPP_
14
15 namespace Bembel {
16 /**
17 * \ingroup LocalOperator
18 * \brief Helper struct which mimics the DiscreteOperatorComputer struct to
19 * assemble sparse matrices of local operators
20 */
21 template <typename Derived>
22 struct DiscreteLocalOperatorComputer {
23 DiscreteLocalOperatorComputer() {}
24 2 static void compute(
25 Eigen::SparseMatrix<typename LinearOperatorTraits<Derived>::Scalar>
26 *disc_op,
27 const Derived &lin_op, const AnsatzSpace<Derived> &ansatz_space) {
28 // Extract numbers
29 2 const auto &super_space = ansatz_space.get_superspace();
30 2 const auto &element_tree = super_space.get_mesh().get_element_tree();
31 2 const auto &number_of_elements = element_tree.get_number_of_elements();
32 2 const auto vector_dimension =
33 getFunctionSpaceVectorDimension<LinearOperatorTraits<Derived>::Form>();
34 2 const auto polynomial_degree = super_space.get_polynomial_degree();
35 2 const auto polynomial_degree_plus_one_squared =
36 2 (polynomial_degree + 1) * (polynomial_degree + 1);
37
38 // Quadrature
39
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 GaussSquare<Constants::maximum_quadrature_degree> GS;
40 2 auto ffield_deg = lin_op.get_FarfieldQuadratureDegree(polynomial_degree);
41
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 auto Q = GS[ffield_deg];
42
43 // Triplets
44 typedef Eigen::Triplet<typename LinearOperatorTraits<Derived>::Scalar> T;
45 2 std::vector<T> tripletList;
46 2 tripletList.reserve(polynomial_degree_plus_one_squared *
47
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 polynomial_degree_plus_one_squared *
48 number_of_elements * vector_dimension);
49
50 // Iterate over elements
51
2/2
✓ Branch 4 taken 12 times.
✓ Branch 5 taken 2 times.
 26 for (auto element = element_tree.cpbegin(); element != element_tree.cpend();
52 12 ++element) {
53 Eigen::Matrix<typename LinearOperatorTraits<Derived>::Scalar,
54 Eigen::Dynamic, Eigen::Dynamic>
55 ✗ intval(vector_dimension * polynomial_degree_plus_one_squared,
56
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 vector_dimension * polynomial_degree_plus_one_squared);
57
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 intval.setZero();
58 // Iterate over quadrature points
59
2/2
✓ Branch 1 taken 108 times.
✓ Branch 2 taken 12 times.
 120 for (auto i = 0; i < Q.w_.size(); ++i) {
60
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 108 SurfacePoint qp;
61
4/8
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 108 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 108 times.
✗ Branch 12 not taken.
 108 super_space.map2surface(*element, Q.xi_.col(i), Q.w_(i), &qp);
62
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 108 lin_op.evaluateIntegrand(super_space, qp, qp, &intval);
63 }
64 // Transform local element matrices to triplets
65
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 12 times.
 24 for (auto i = 0; i < vector_dimension; ++i)
66
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 12 times.
 24 for (auto j = 0; j < vector_dimension; ++j)
67
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 12 times.
 60 for (auto ii = 0; ii < polynomial_degree_plus_one_squared; ++ii)
68
2/2
✓ Branch 0 taken 192 times.
✓ Branch 1 taken 48 times.
 240 for (auto jj = 0; jj < polynomial_degree_plus_one_squared; ++jj)
69
1/2
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
 192 tripletList.push_back(
70 ✗ T(polynomial_degree_plus_one_squared *
71 192 (j * number_of_elements + element->id_) +
72 jj,
73 192 polynomial_degree_plus_one_squared *
74 384 (i * number_of_elements + element->id_) +
75 ii,
76 192 intval(j * polynomial_degree_plus_one_squared + jj,
77
1/2
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
 192 i * polynomial_degree_plus_one_squared + ii)));
78 }
79 2 disc_op->resize(vector_dimension * polynomial_degree_plus_one_squared *
80 number_of_elements,
81
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 vector_dimension * polynomial_degree_plus_one_squared *
82 number_of_elements);
83
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 disc_op->setFromTriplets(tripletList.begin(), tripletList.end());
84 2 const auto &projector = ansatz_space.get_transformation_matrix();
85
4/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 2 disc_op[0] = projector.transpose() * (disc_op[0] * projector);
86 4 return;
87 2 }
88 };
89
90 /**
91 * \ingroup LocalOperator
92 * \brief DiscreteLocalOperator
93 * Specialization of the DiscreteOperator class for sparse matrices
94 * and an intgrator for local operators
95 **/
96 template <typename Derived>
97 using DiscreteLocalOperator = DiscreteOperator<
98 Eigen::SparseMatrix<typename LinearOperatorTraits<Derived>::Scalar>,
99 Derived>;
100
101 /**
102 * \brief Helper struct that is used in order to partially specialise the
103 * compute routine of DiscreteOperator for the Eigen::H2Matrix format
104 */
105 template <typename Derived>
106 struct DiscreteOperatorComputer<
107 Eigen::SparseMatrix<typename LinearOperatorTraits<Derived>::Scalar>,
108 Derived> {
109 DiscreteOperatorComputer() {}
110 2 static void compute(
111 Eigen::SparseMatrix<typename LinearOperatorTraits<Derived>::Scalar>
112 *disc_op,
113 const Derived &lin_op, const AnsatzSpace<Derived> &ansatz_space) {
114 2 DiscreteLocalOperatorComputer<Derived>::compute(disc_op, lin_op,
115 ansatz_space);
116 2 return;
117 }
118 };
119
120 } // namespace Bembel
121 #endif // BEMBEL_SRC_LINEAROPERATOR_DISCRETELOCALOPERATOR_HPP_
122