GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/AnsatzSpace/Projector.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 154 156 98.7%
Functions: 77 77 100.0%
Branches: 106 164 64.6%
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_ANSATZSPACE_PROJECTOR_HPP_
12 #define BEMBEL_SRC_ANSATZSPACE_PROJECTOR_HPP_
13
14 namespace Bembel {
15 namespace ProjectorRoutines {
16 /**
17 * \brief Helper struct for assembling the Projector for the different cases.
18 */
19 template <typename Derived, unsigned int DF>
20 struct projector_matrixmaker_ {
21 static Eigen::SparseMatrix<double> makeMatrix(
22 const Bembel::SuperSpace<Derived>& super_space,
23 const int knotrepetition) {
24 assert(false && "This needs to be specialized");
25 return Eigen::SparseMatrix<double>(1, 1);
26 }
27 };
28 } // namespace ProjectorRoutines
29
30 /**
31 * \ingroup AnsatzSpace
32 * \brief The Projector provides routines to assemble smooth B-Splines on each
33 * patch.
34 *
35 * This class solves local interpolation problems on each patch to assemble
36 * B-splines from the Berstein basis. The linear combination of the Bernstein
37 * basis functions is stored in a sparse transformation matrix.
38 */
39 template <typename Derived>
40 class Projector {
41 public:
42 //////////////////////////////////////////////////////////////////////////////
43 /// constructors
44 //////////////////////////////////////////////////////////////////////////////
45 /**
46 * \brief Default constructor
47 */
48 Projector() {}
49 /**
50 * \brief Constructor for the Projector.
51 *
52 * This constructor initializes a Projector object with the provided
53 * parameters.
54 *
55 * \param super_space The SuperSpace to handle basis functions.
56 * \param knot_repetition The number of repetitions of knots in the space.
57 */
58 344 Projector(const SuperSpace<Derived>& super_space, const int knot_repetition) {
59
1/2
✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
 344 init_Projector(super_space, knot_repetition);
60 344 }
61 //////////////////////////////////////////////////////////////////////////////
62 /// init
63 //////////////////////////////////////////////////////////////////////////////
64 /**
65 * \brief Initializes the Projector.
66 *
67 * This function initializes the member variables of the Projector
68 *
69 * \param super_space The SuperSpace to handle basis functions.
70 * \param knot_repetition The number of repetitions of knots in the space.
71 */
72 344 void init_Projector(const SuperSpace<Derived>& super_space,
73 const int knot_repetition) {
74 344 knot_repetition_ = knot_repetition;
75
1/2
✓ Branch 2 taken 180 times.
✗ Branch 3 not taken.
 344 projector_ = makeProjectionMatrix(super_space, knot_repetition);
76 344 dofs_before_projector_ = projector_.rows();
77 344 dofs_after_projector_ = projector_.cols();
78 344 return;
79 }
80 //////////////////////////////////////////////////////////////////////////////
81 /// getter
82 //////////////////////////////////////////////////////////////////////////////
83 /**
84 * \brief Return knot repetition.
85 *
86 * \return The number of repetitions of knots in the space.
87 */
88 392 int get_knot_repetition() const { return knot_repetition_; }
89 /**
90 * \brief This function returns the transformation matrix to assemble smooth
91 * tensor product B-splines.
92 *
93 * In the assembly process of this Projection matrix an local interpolation
94 * problem is solved on each element. This matrix needs to be applied to the
95 * Bernstein basis to assemble smooth tensor product B-splines.
96 *
97 * \return The Projector matrix a tall thin Sparse transformation matrix.
98 */
99 344 const Eigen::SparseMatrix<double>& get_projection_matrix() {
100 344 return projector_;
101 }
102 /**
103 * \brief Return number of smooth B-splines.
104 *
105 * \return Number of DOFs after applying the transformation matrix.
106 */
107 533 int get_dofs_after_projector() const { return dofs_after_projector_; }
108 //////////////////////////////////////////////////////////////////////////////
109 /// private members
110 //////////////////////////////////////////////////////////////////////////////
111 private:
112 344 Eigen::SparseMatrix<double> makeProjectionMatrix(
113 const SuperSpace<Derived>& super_space, const int knotrepetition) {
114 return ProjectorRoutines::projector_matrixmaker_<
115 Derived,
116 LinearOperatorTraits<Derived>::Form>::makeMatrix(super_space,
117 344 knotrepetition);
118 }
119 // we declare functionality which has not been implemented (yet)
120 // to be private
121 int knot_repetition_;
122 int dofs_before_projector_;
123 int dofs_after_projector_;
124 Eigen::SparseMatrix<double> projector_;
125 Projector(const Projector<Derived>& other);
126 Projector(Projector<Derived>&& other);
127 Projector& operator=(const Projector<Derived>& other);
128 Projector& operator=(Projector<Derived>&& other);
129 //////////////////////////////////////////////////////////////////////////////
130 };
131
132 namespace ProjectorRoutines {
133
134 /**
135 * \brief Helper struct
136 *
137 * \todo Use Eigen::Triplets.
138 */
139 struct _proj_info {
140 std::vector<double> vals;
141 std::vector<int> rows;
142 std::vector<int> cols;
143 int num_c_dof;
144 int num_dc_dof;
145 };
146
147 template <typename Derived>
148 /**
149 * \brief This function solves on each element the interpolation problem to
150 * assemble smooth tensor product B-splines.
151 *
152 * This function returns the entries of the matrix which transforms Bernstein
153 * polynomials to tensor product B-splines defined by the polynomial degree in x
154 * and y direction. The local interpolation mask is filled with the DeBoor
155 * algorithm and the system is solved with the Eigen::FullPivLU decomposition.
156 *
157 * \param super_space Reference to the SuperSpace handling basis functions.
158 * \param pp1x Polynomial degree in x direction plus 1.
159 * \param pp1y Polynomial degree in y direction plus 1.
160 * \param knotrepetition_in Knot repetition of the knot vector.
161 *
162 * \return Struct Number of DOFs and indices of the matrix entries and the
163 * coefficients of the linear combination assembling smooth B-spline.
164 */
165 499 inline _proj_info makeLocalProjectionTriplets(
166 const SuperSpace<Derived>& super_space, const int pp1x, const int pp1y,
167 const int knotrepetition_in) {
168 499 const int M = super_space.get_refinement_level();
169 499 const int maximal_polynomial_degree = std::max(pp1x, pp1y);
170
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 assert(maximal_polynomial_degree == super_space.get_polynomial_degree() + 1 &&
171 "superSpace not suitable for desired ansatz space -> polynomial "
172 "degree mismatch");
173 499 const int minp = std::min(pp1x, pp1y);
174 499 const int knotrepetition = (maximal_polynomial_degree <= knotrepetition_in)
175
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 237 times.
 499 ? minp
176 : knotrepetition_in;
177 499 const int n = (1 << M);
178 499 const int patch_number = super_space.get_number_of_patches();
179
180 // At the moment, we allow a difference of one only, i.e., we start our de
181 // Rham sequence with a space of the same degree in every TP-direction.
182
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 259 times.
 499 assert(std::abs(maximal_polynomial_degree - minp) <= 1);
183
184 // Now, we construct the continuous spaces which will be the preimage of the
185 // projector. n-1 is passed, since n = number_elements but n-1 = number of
186 // interior knots.
187 499 std::vector<double> c_space_knot_x =
188 Spl::MakeUniformKnotVector(pp1x, n - 1, knotrepetition);
189 499 std::vector<double> c_space_knot_y =
190 Spl::MakeUniformKnotVector(pp1y, n - 1, knotrepetition);
191 499 const int c_space_dim_x = c_space_knot_x.size() - pp1x;
192 499 const int c_space_dim_y = c_space_knot_y.size() - pp1y;
193 499 const int c_space_dim = c_space_dim_x * c_space_dim_y;
194
195 499 _proj_info out;
196
197 499 out.cols.reserve(c_space_dim * maximal_polynomial_degree *
198
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 maximal_polynomial_degree * patch_number);
199 499 out.vals.reserve(c_space_dim * maximal_polynomial_degree *
200
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 maximal_polynomial_degree * patch_number);
201 499 out.rows.reserve(c_space_dim * maximal_polynomial_degree *
202
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 maximal_polynomial_degree * patch_number);
203
204 499 std::vector<double> mask =
205 Spl::MakeInterpolationMask(maximal_polynomial_degree);
206 499 const int masksize = mask.size();
207
208 // Now we assemble our system which we use for the interpolation
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 259 times.
 499 assert(masksize == maximal_polynomial_degree &&
210 "projector.cpp: System needs to be square");
211 499 Eigen::Matrix<double, -1, -1> system(
212 ✗ masksize * masksize,
213
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 maximal_polynomial_degree * maximal_polynomial_degree);
214
215 // Here, we suddenly use the degree for basisevaluation, i.e.,
216 // maximal_polynomial_degree-1. This is confusing, but correct and tested.
217 {
218 double vals_y[Constants::MaxP + 1];
219 double vals_x[Constants::MaxP + 1];
220
2/2
✓ Branch 0 taken 868 times.
✓ Branch 1 taken 259 times.
 2207 for (int iy = 0; iy < maximal_polynomial_degree; ++iy) {
221
1/2
✓ Branch 2 taken 868 times.
✗ Branch 3 not taken.
 1708 Bembel::Basis::ShapeFunctionHandler::evalBasis(
222 maximal_polynomial_degree - 1, vals_y, mask[iy]);
223
2/2
✓ Branch 0 taken 3286 times.
✓ Branch 1 taken 868 times.
 8234 for (int ix = 0; ix < maximal_polynomial_degree; ++ix) {
224
1/2
✓ Branch 2 taken 3286 times.
✗ Branch 3 not taken.
 6526 Bembel::Basis::ShapeFunctionHandler::evalBasis(
225 maximal_polynomial_degree - 1, vals_x, mask[ix]);
226
2/2
✓ Branch 0 taken 13522 times.
✓ Branch 1 taken 3286 times.
 33488 for (int jy = 0; jy < maximal_polynomial_degree; ++jy) {
227
2/2
✓ Branch 0 taken 58834 times.
✓ Branch 1 taken 13522 times.
 144476 for (int jx = 0; jx < maximal_polynomial_degree; ++jx) {
228 235028 system(iy * masksize + ix, jy * maximal_polynomial_degree + jx) =
229
1/2
✓ Branch 1 taken 58834 times.
✗ Branch 2 not taken.
 117514 vals_x[jx] * vals_y[jy];
230 }
231 }
232 }
233 }
234 }
235
236 // Do LU-Decomposition of
237 // systemsuper_space.get_mesh().get_element_tree().begin()
238
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 const Eigen::PartialPivLU<Eigen::Matrix<double, -1, -1>> pplu(system);
239
240 499 unsigned int element_number = 0;
241
242 46839 for (auto element = super_space.get_mesh().get_element_tree().cpbegin();
243
2/2
✓ Branch 4 taken 25604 times.
✓ Branch 5 taken 259 times.
 46839 element != super_space.get_mesh().get_element_tree().cpend();
244 46340 ++element) {
245 // s = x
246
1/2
✓ Branch 2 taken 25604 times.
✗ Branch 3 not taken.
 46340 const double pos_x = element->llc_(0);
247
1/2
✓ Branch 2 taken 25604 times.
✗ Branch 3 not taken.
 46340 const double pos_y = element->llc_(1);
248
249 46340 const double h = element->get_h();
250 46340 const double mid_x = pos_x + (.5 * h);
251 46340 const double mid_y = pos_y + (.5 * h);
252
253 // The following block of code computes the indeces of all continuous basis
254 // functions which have a support on the element of index 'element'
255 46340 std::vector<double> nonzero_dofs_x;
256 46340 std::vector<double> nonzero_dofs_y;
257
1/2
✓ Branch 1 taken 25604 times.
✗ Branch 2 not taken.
 46340 nonzero_dofs_x.reserve(pp1x);
258
1/2
✓ Branch 1 taken 25604 times.
✗ Branch 2 not taken.
 46340 nonzero_dofs_y.reserve(pp1y);
259
2/2
✓ Branch 0 taken 452867 times.
✓ Branch 1 taken 25604 times.
 686251 for (int dof_y = 0; dof_y < c_space_dim_y; ++dof_y) {
260
4/4
✓ Branch 1 taken 261819 times.
✓ Branch 2 taken 191048 times.
✓ Branch 3 taken 70771 times.
✓ Branch 4 taken 382096 times.
 1028098 if ((c_space_knot_y[dof_y] <= mid_y) &&
261
2/2
✓ Branch 1 taken 70771 times.
✓ Branch 2 taken 191048 times.
 388187 (c_space_knot_y[dof_y + pp1y] >= mid_y)) {
262
1/2
✓ Branch 1 taken 70771 times.
✗ Branch 2 not taken.
 136463 nonzero_dofs_y.push_back(dof_y);
263 }
264 }
265
2/2
✓ Branch 0 taken 452867 times.
✓ Branch 1 taken 25604 times.
 686251 for (int dof_x = 0; dof_x < c_space_dim_x; ++dof_x) {
266
4/4
✓ Branch 1 taken 261819 times.
✓ Branch 2 taken 191048 times.
✓ Branch 3 taken 70771 times.
✓ Branch 4 taken 382096 times.
 1028098 if ((c_space_knot_x[dof_x] <= mid_x) &&
267
2/2
✓ Branch 1 taken 70771 times.
✓ Branch 2 taken 191048 times.
 388187 (c_space_knot_x[dof_x + pp1x] >= mid_x)) {
268
1/2
✓ Branch 1 taken 70771 times.
✗ Branch 2 not taken.
 136463 nonzero_dofs_x.push_back(dof_x);
269 }
270 }
271
1/2
✓ Branch 1 taken 25604 times.
✗ Branch 2 not taken.
 46340 nonzero_dofs_x.shrink_to_fit();
272
1/2
✓ Branch 1 taken 25604 times.
✗ Branch 2 not taken.
 46340 nonzero_dofs_y.shrink_to_fit();
273
274 // We now loop over all basis functions wich are non-zero on the element
275
2/2
✓ Branch 5 taken 70771 times.
✓ Branch 6 taken 25604 times.
 182803 for (auto dof_y : nonzero_dofs_y) {
276
2/2
✓ Branch 13 taken 237146 times.
✓ Branch 14 taken 70771 times.
 605357 for (auto dof_x : nonzero_dofs_x) {
277 // The next few lines are to check if the element is within the support
278 // of the basis function given by dof_x and dof_y.
279
1/2
✓ Branch 2 taken 237146 times.
✗ Branch 3 not taken.
 468894 std::vector<double> local_mask_x(masksize);
280
1/2
✓ Branch 2 taken 237146 times.
✗ Branch 3 not taken.
 468894 std::vector<double> local_mask_y(masksize);
281
2/2
✓ Branch 0 taken 981590 times.
✓ Branch 1 taken 237146 times.
 2425832 for (int i = 0; i < masksize; ++i) {
282 1956938 local_mask_x[i] = pos_x + (1.0 / n) * mask[i];
283 1956938 local_mask_y[i] = pos_y + (1.0 / n) * mask[i];
284
2/4
✓ Branch 1 taken 981590 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 981590 times.
✗ Branch 5 not taken.
 1956938 assert(local_mask_x[i] < 1 && local_mask_x[i] > 0);
285
2/4
✓ Branch 1 taken 981590 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 981590 times.
✗ Branch 5 not taken.
 1956938 assert(local_mask_y[i] < 1 && local_mask_y[i] > 0);
286 }
287
288 // build unit vectors
289
1/2
✓ Branch 2 taken 237146 times.
✗ Branch 3 not taken.
 468894 std::vector<double> c_coefs_x(c_space_dim_x, 0.0);
290
1/2
✓ Branch 2 taken 237146 times.
✗ Branch 3 not taken.
 468894 std::vector<double> c_coefs_y(c_space_dim_y, 0.0);
291 468894 c_coefs_x[dof_x] = 1.;
292 468894 c_coefs_y[dof_y] = 1.;
293
294 // evaluate rhs for interpolation problem
295 468894 std::vector<double> vals_x =
296 Spl::DeBoor(c_coefs_x, c_space_knot_x, local_mask_x);
297 468894 std::vector<double> vals_y =
298 Spl::DeBoor(c_coefs_y, c_space_knot_y, local_mask_y);
299
300 // assemble and solve interpolation problem
301
1/2
✓ Branch 1 taken 237146 times.
✗ Branch 2 not taken.
 468894 Eigen::VectorXd rhs(masksize * masksize);
302
2/2
✓ Branch 0 taken 981590 times.
✓ Branch 1 taken 237146 times.
 2425832 for (int iy = 0; iy < masksize; ++iy)
303
2/2
✓ Branch 0 taken 4306874 times.
✓ Branch 1 taken 981590 times.
 10562756 for (int ix = 0; ix < masksize; ++ix)
304
1/2
✓ Branch 3 taken 4306874 times.
✗ Branch 4 not taken.
 8605818 rhs(iy * masksize + ix) = vals_x[ix] * vals_y[iy];
305
2/4
✓ Branch 1 taken 237146 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 237146 times.
✗ Branch 5 not taken.
 468894 Eigen::VectorXd sol = pplu.solve(rhs);
306
307 // insert entries into helperstruct
308 9074712 for (int i = 0;
309
2/2
✓ Branch 0 taken 4306874 times.
✓ Branch 1 taken 237146 times.
 9074712 i < maximal_polynomial_degree * maximal_polynomial_degree; ++i) {
310
3/4
✓ Branch 1 taken 4306874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1552022 times.
✓ Branch 5 taken 2754852 times.
 8605818 if (std::abs(sol(i)) > Constants::projector_tolerance) {
311 3098234 out.cols.push_back(element->patch_ * (c_space_dim) +
312
1/2
✓ Branch 1 taken 1552022 times.
✗ Branch 2 not taken.
 3098234 (c_space_dim_x * dof_y) + dof_x);
313 3098234 out.rows.push_back((element_number * maximal_polynomial_degree *
314 ✗ maximal_polynomial_degree) +
315
1/2
✓ Branch 1 taken 1552022 times.
✗ Branch 2 not taken.
 3098234 i);
316
2/4
✓ Branch 1 taken 1552022 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1552022 times.
✗ Branch 5 not taken.
 3098234 out.vals.push_back(sol(i));
317 }
318 }
319 }
320 }
321 46340 ++element_number;
322 }
323
324 499 out.num_c_dof = c_space_dim_x * c_space_dim_y * patch_number;
325 499 out.num_dc_dof = maximal_polynomial_degree * maximal_polynomial_degree * n *
326 499 n * patch_number;
327
328
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 out.cols.shrink_to_fit();
329
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 out.rows.shrink_to_fit();
330
1/2
✓ Branch 1 taken 259 times.
✗ Branch 2 not taken.
 499 out.vals.shrink_to_fit();
331
332
2/4
✓ Branch 2 taken 259 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 259 times.
✗ Branch 7 not taken.
 499 assert((out.vals.size() == out.rows.size()) &&
333 (out.rows.size() == out.cols.size()) &&
334 "projector.cpp: you made a mistake, try again.");
335
336 998 return out;
337 499 }
338
339 /**
340 * \brief Specification of the struct for the scalar continuous case.
341 *
342 * Do the same as in the discontinuous case but with an assert p >= 1.
343 */
344 template <typename Derived>
345 struct projector_matrixmaker_<Derived, DifferentialForm::Continuous> {
346 66 static Eigen::SparseMatrix<double> makeMatrix(
347 const SuperSpace<Derived>& super_space, const int knotrepetition) {
348 66 const int P = super_space.get_polynomial_degree();
349
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 66 times.
 66 assert(P > 0 && "P must be 1 or larger for this type of discrete space");
350
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 66 times.
 66 assert(knotrepetition <= P &&
351 "Knot repetition must be smaller than P for this type of discrete "
352 "space");
353 // The matrices coincide with the 2 case, up to the limitations above
354 66 Eigen::SparseMatrix<double> local_matrix = projector_matrixmaker_<
355 Derived, DifferentialForm::Discontinuous>::makeMatrix(super_space,
356 knotrepetition);
357
358 66 return local_matrix;
359 }
360 };
361
362 /**
363 * \brief Specification of the struct for the vector Div conforming case.
364 *
365 * The Maxwell case: This function assembles a vector valued smooth tensor
366 * product B-spline basis. Hereby the first component has the polynomial degree
367 * p in x direction and p - 1 in y direction. For the second component the
368 * polynomial degree for x and y direction are switched.
369 */
370 template <typename Derived>
371 struct projector_matrixmaker_<Derived, DifferentialForm::DivConforming> {
372 79 static Eigen::SparseMatrix<double> makeMatrix(
373 const SuperSpace<Derived>& super_space, const int knotrepetition) {
374
1/2
✓ Branch 1 taken 79 times.
✗ Branch 2 not taken.
 79 assert(super_space.get_polynomial_degree() > 0 &&
375 "P must be 1 or larger for this type of discrete space");
376
1/2
✓ Branch 1 taken 79 times.
✗ Branch 2 not taken.
 79 assert(knotrepetition <= super_space.get_polynomial_degree() &&
377 "Knot repetition must be smaller than P for this type of discrete "
378 "space");
379
1/2
✓ Branch 2 taken 79 times.
✗ Branch 3 not taken.
 158 const auto info1 = ProjectorRoutines::makeLocalProjectionTriplets<Derived>(
380 79 super_space, super_space.get_polynomial_degree() + 1,
381 super_space.get_polynomial_degree(), knotrepetition);
382
1/2
✓ Branch 2 taken 79 times.
✗ Branch 3 not taken.
 79 const auto info2 = ProjectorRoutines::makeLocalProjectionTriplets<Derived>(
383 super_space, super_space.get_polynomial_degree(),
384 79 super_space.get_polynomial_degree() + 1, knotrepetition);
385 79 const int size1 = info1.vals.size();
386 79 const int size2 = info2.vals.size();
387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 79 times.
 79 assert(size1 == size2);
388 79 std::vector<Eigen::Triplet<double>> trips;
389
2/2
✓ Branch 0 taken 489744 times.
✓ Branch 1 taken 79 times.
 489823 for (int k = 0; k < size1; ++k) {
390
1/2
✓ Branch 1 taken 489744 times.
✗ Branch 2 not taken.
 489744 trips.push_back(
391 979488 Eigen::Triplet<double>(info1.rows[k], info1.cols[k], info1.vals[k]));
392 }
393
2/2
✓ Branch 0 taken 489744 times.
✓ Branch 1 taken 79 times.
 489823 for (int k = 0; k < size2; ++k) {
394
1/2
✓ Branch 2 taken 489744 times.
✗ Branch 3 not taken.
 489744 trips.push_back(Eigen::Triplet<double>(info2.rows[k] + info1.num_dc_dof,
395 979488 info2.cols[k] + info1.num_c_dof,
396 489744 info2.vals[k]));
397 }
398 79 Eigen::SparseMatrix<double> local_matrix(
399
1/2
✓ Branch 1 taken 79 times.
✗ Branch 2 not taken.
 79 info1.num_dc_dof + info2.num_dc_dof, info1.num_c_dof + info2.num_c_dof);
400
1/2
✓ Branch 3 taken 79 times.
✗ Branch 4 not taken.
 79 local_matrix.setFromTriplets(trips.begin(), trips.end());
401
402 158 return local_matrix;
403 79 }
404 };
405
406 /**
407 * \brief Specification of the struct for the discontinuous case.
408 *
409 * This functions assembles scalar smooth tensor product B-splines according to
410 * the knot repetition. The polynomial degree in x and y direction is the same.
411 */
412 template <typename Derived>
413 struct projector_matrixmaker_<Derived, DifferentialForm::Discontinuous> {
414 101 static Eigen::SparseMatrix<double> makeMatrix(
415 const SuperSpace<Derived>& super_space, const int knotrepetition) {
416
1/2
✓ Branch 1 taken 101 times.
✗ Branch 2 not taken.
 101 const auto info = ProjectorRoutines::makeLocalProjectionTriplets<Derived>(
417 101 super_space, super_space.get_polynomial_degree() + 1,
418 101 super_space.get_polynomial_degree() + 1, knotrepetition);
419 101 const int size = info.vals.size();
420 101 std::vector<Eigen::Triplet<double>> trips;
421
2/2
✓ Branch 0 taken 572534 times.
✓ Branch 1 taken 101 times.
 572635 for (int k = 0; k < size; ++k) {
422
1/2
✓ Branch 1 taken 572534 times.
✗ Branch 2 not taken.
 572534 trips.push_back(
423 1145068 Eigen::Triplet<double>(info.rows[k], info.cols[k], info.vals[k]));
424 }
425
1/2
✓ Branch 1 taken 101 times.
✗ Branch 2 not taken.
 101 Eigen::SparseMatrix<double> local_matrix(info.num_dc_dof, info.num_c_dof);
426
1/2
✓ Branch 3 taken 101 times.
✗ Branch 4 not taken.
 101 local_matrix.setFromTriplets(trips.begin(), trips.end());
427
428 202 return local_matrix;
429 101 }
430 };
431 } // namespace ProjectorRoutines
432
433 } // namespace Bembel
434
435 #endif // BEMBEL_SRC_ANSATZSPACE_PROJECTOR_HPP_
436