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 | #ifndef BEMBEL_SRC_LINEAROPERATOR_DISCRETEOPERATOR_HPP_ | ||
12 | #define BEMBEL_SRC_LINEAROPERATOR_DISCRETEOPERATOR_HPP_ | ||
13 | |||
14 | namespace Bembel { | ||
15 | /** | ||
16 | * \ingroup LinearOperator | ||
17 | * \brief Helper struct that is used in order to partially specialise the | ||
18 | * compute routine of DiscreteOperator for different types of | ||
19 | * output formats | ||
20 | */ | ||
21 | template <typename MatrixFormat, typename Derived> | ||
22 | struct DiscreteOperatorComputer {}; | ||
23 | /** | ||
24 | * \brief Helper struct that is used in order to partially specialise the | ||
25 | * compute routine of DiscreteOperator for the Eigen::MatrixXd format | ||
26 | */ | ||
27 | template <typename Derived, typename Scalar> | ||
28 | struct DiscreteOperatorComputer< | ||
29 | Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>, Derived> { | ||
30 | DiscreteOperatorComputer() {} | ||
31 | 3 | static void compute( | |
32 | Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> *disc_op, | ||
33 | const Derived &lin_op, const AnsatzSpace<Derived> &ansatz_space) { | ||
34 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | GaussSquare<Constants::maximum_quadrature_degree> GS; |
35 | 3 | const SuperSpace<Derived> &super_space = ansatz_space.get_superspace(); | |
36 | 3 | const ElementTree &element_tree = super_space.get_mesh().get_element_tree(); | |
37 | 3 | auto number_of_elements = element_tree.get_number_of_elements(); | |
38 | 3 | const auto vector_dimension = | |
39 | getFunctionSpaceVectorDimension<LinearOperatorTraits<Derived>::Form>(); | ||
40 | 3 | auto polynomial_degree = super_space.get_polynomial_degree(); | |
41 | 3 | auto polynomial_degree_plus_one_squared = | |
42 | 3 | (polynomial_degree + 1) * (polynomial_degree + 1); | |
43 | 3 | auto ffield_deg = lin_op.get_FarfieldQuadratureDegree(polynomial_degree); | |
44 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
3 | auto ffield_qnodes = |
45 | DuffyTrick::computeFfieldQnodes(super_space, GS[ffield_deg]); | ||
46 | 3 | disc_op->resize(vector_dimension * polynomial_degree_plus_one_squared * | |
47 | number_of_elements, | ||
48 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | vector_dimension * polynomial_degree_plus_one_squared * |
49 | number_of_elements); | ||
50 | #pragma omp parallel | ||
51 | { | ||
52 | #pragma omp single | ||
53 | { | ||
54 | 3 | for (auto element1 = element_tree.cpbegin(); | |
55 |
2/2✓ Branch 3 taken 126 times.
✓ Branch 4 taken 3 times.
|
129 | element1 != element_tree.cpend(); ++element1) |
56 | 126 | for (auto element2 = element_tree.cpbegin(); | |
57 |
2/2✓ Branch 3 taken 9828 times.
✓ Branch 4 taken 126 times.
|
9954 | element2 != element_tree.cpend(); ++element2) { |
58 | #pragma omp task firstprivate(element1, element2) | ||
59 | { | ||
60 | 9828 | Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> intval( | |
61 | ✗ | vector_dimension * polynomial_degree_plus_one_squared, | |
62 |
1/2✓ Branch 1 taken 9828 times.
✗ Branch 2 not taken.
|
9828 | vector_dimension * polynomial_degree_plus_one_squared); |
63 |
1/2✓ Branch 1 taken 9828 times.
✗ Branch 2 not taken.
|
9828 | DuffyTrick::evaluateBilinearForm( |
64 | 9828 | lin_op, super_space, *element1, *element2, GS, | |
65 | 9828 | ffield_qnodes[element1->id_], ffield_qnodes[element2->id_], | |
66 | &intval); | ||
67 |
2/2✓ Branch 0 taken 9864 times.
✓ Branch 1 taken 9828 times.
|
19692 | for (auto i = 0; i < vector_dimension; ++i) |
68 |
2/2✓ Branch 0 taken 9936 times.
✓ Branch 1 taken 9864 times.
|
19800 | for (auto j = 0; j < vector_dimension; ++j) |
69 | ✗ | disc_op->block(polynomial_degree_plus_one_squared * | |
70 |
1/2✓ Branch 1 taken 9936 times.
✗ Branch 2 not taken.
|
9936 | (i * number_of_elements + element1->id_), |
71 | 9936 | polynomial_degree_plus_one_squared * | |
72 | 9936 | (j * number_of_elements + element2->id_), | |
73 | polynomial_degree_plus_one_squared, | ||
74 |
1/2✓ Branch 1 taken 9936 times.
✗ Branch 2 not taken.
|
9936 | polynomial_degree_plus_one_squared) = |
75 | 9936 | intval.block(i * polynomial_degree_plus_one_squared, | |
76 |
1/2✓ Branch 1 taken 9936 times.
✗ Branch 2 not taken.
|
9936 | j * polynomial_degree_plus_one_squared, |
77 | polynomial_degree_plus_one_squared, | ||
78 | polynomial_degree_plus_one_squared); | ||
79 | 9828 | } | |
80 | } | ||
81 | } | ||
82 | } | ||
83 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | auto projector = ansatz_space.get_transformation_matrix(); |
84 |
4/8✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
|
3 | disc_op[0] = projector.transpose() * (disc_op[0] * projector); |
85 | 6 | return; | |
86 | 3 | } | |
87 | }; | ||
88 | /** | ||
89 | * \brief Helper struct that is used in order to partially specialise the | ||
90 | * compute routine of DiscreteOperator for the Eigen::H2Matrix format | ||
91 | */ | ||
92 | template <typename Derived> | ||
93 | struct DiscreteOperatorComputer< | ||
94 | Eigen::H2Matrix<typename LinearOperatorTraits<Derived>::Scalar>, Derived> { | ||
95 | DiscreteOperatorComputer() {} | ||
96 | 5 | static void compute( | |
97 | Eigen::H2Matrix<typename LinearOperatorTraits<Derived>::Scalar> *disc_op, | ||
98 | const Derived &lin_op, const AnsatzSpace<Derived> &ansatz_space) { | ||
99 | 5 | disc_op->init_H2Matrix(lin_op, ansatz_space); | |
100 | 5 | return; | |
101 | } | ||
102 | }; | ||
103 | /** | ||
104 | * \brief DiscreteOperator | ||
105 | * \todo Add a documentation | ||
106 | */ | ||
107 | template <typename MatrixFormat, typename Derived> | ||
108 | class DiscreteOperator { | ||
109 | public: | ||
110 | ////////////////////////////////////////////////////////////////////////////// | ||
111 | // constructors | ||
112 | ////////////////////////////////////////////////////////////////////////////// | ||
113 | DiscreteOperator() {} | ||
114 |
3/5✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
|
12 | explicit DiscreteOperator(const AnsatzSpace<Derived> &ansatz_space) { |
115 |
1/2✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
|
12 | init_DiscreteOperator(ansatz_space); |
116 | 12 | } | |
117 | ////////////////////////////////////////////////////////////////////////////// | ||
118 | // init_DiscreteOperator | ||
119 | ////////////////////////////////////////////////////////////////////////////// | ||
120 | 12 | void init_DiscreteOperator(const AnsatzSpace<Derived> &ansatz_space) { | |
121 |
1/2✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
12 | ansatz_space_ = ansatz_space; |
122 | 12 | return; | |
123 | } | ||
124 | ////////////////////////////////////////////////////////////////////////////// | ||
125 | // compute | ||
126 | ////////////////////////////////////////////////////////////////////////////// | ||
127 | 12 | inline void compute() { | |
128 | 12 | DiscreteOperatorComputer<MatrixFormat, Derived>::compute(&disc_op_, lin_op_, | |
129 | 12 | ansatz_space_); | |
130 | 12 | return; | |
131 | } | ||
132 | ////////////////////////////////////////////////////////////////////////////// | ||
133 | // getter | ||
134 | ////////////////////////////////////////////////////////////////////////////// | ||
135 | 13 | const MatrixFormat &get_discrete_operator() const { return disc_op_; } | |
136 | const Derived &get_linear_operator() const { return lin_op_; } | ||
137 | 3 | Derived &get_linear_operator() { return lin_op_; } | |
138 | ////////////////////////////////////////////////////////////////////////////// | ||
139 | // private member variables | ||
140 | ////////////////////////////////////////////////////////////////////////////// | ||
141 | private: | ||
142 | Derived lin_op_; | ||
143 | MatrixFormat disc_op_; | ||
144 | AnsatzSpace<Derived> ansatz_space_; | ||
145 | }; | ||
146 | |||
147 | } // namespace Bembel | ||
148 | #endif // BEMBEL_SRC_LINEAROPERATOR_DISCRETEOPERATOR_HPP_ | ||
149 |