GCC Code Coverage Report

Directory: Bembel/src/
File: Geometry/Patch.hpp
Date: 2024-12-18 07:36:36
Exec Total Coverage
Lines: 178 178 100.0%
Functions: 9 9 100.0%
Branches: 112 166 67.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 #ifndef BEMBEL_SRC_GEOMETRY_PATCH_HPP_
12 #define BEMBEL_SRC_GEOMETRY_PATCH_HPP_
13 namespace Bembel {
14
15 /**
16 * \ingroup Geometry
17 * \class Patch
18 * \brief handles a single patch
19 */
20 class Patch {
21 public:
22 //////////////////////////////////////////////////////////////////////////////
23 /// constructors
24 //////////////////////////////////////////////////////////////////////////////
25 /**
26 * \brief Default constructor.
27 */
28 98 Patch() {}
29 /**
30 * \brief Constructor of a Patch.
31 *
32 * A patch is defined by a set of control points and knot vectors in the x and
33 * y directions. It is used to create complex surfaces in computer graphics
34 * and CAD applications.
35 *
36 * \param control_points std::vector<Eigen::MatrixXd> {x, y, z, w} where the
37 * rows are the control points in x direction.
38 * \param knots_x The knot vector
39 * in the x direction.
40 * \param knots_y The knot vector in the y direction.
41 */
42 Patch(const std::vector<Eigen::Matrix<double, -1, -1>> &control_points,
43 const std::vector<double> &knots_x,
44 const std::vector<double> &knots_y) {
45 init_Patch(control_points, knots_x, knots_y);
46 }
47 //////////////////////////////////////////////////////////////////////////////
48 /// init
49 //////////////////////////////////////////////////////////////////////////////
50 /**
51 * \brief Initializes a Patch.
52 *
53 * A patch is defined by a set of control points and knot vectors in the x and
54 * y directions. It is used to create complex surfaces in computer graphics
55 * and CAD applications.
56 *
57 * \param xyzw std::vector<Eigen::MatrixXd> {x, y, z, w} where the
58 * rows are the control points in x direction.
59 * \param x_knots The knot vector
60 * in the x direction.
61 * \param y_knots The knot vector in the y direction.
62 */
63 94 inline void init_Patch(const std::vector<Eigen::Matrix<double, -1, -1>> &xyzw,
64 const std::vector<double> &x_knots,
65 const std::vector<double> &y_knots) {
66
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 94 times.
 94 assert(xyzw.size() == 4);
67 94 const int xyzw_cols = xyzw[0].cols();
68 94 const int xyzw_rows = xyzw[0].rows();
69 94 unique_knots_x_ = Spl::ExtractUniqueKnotVector(x_knots);
70 94 unique_knots_y_ = Spl::ExtractUniqueKnotVector(y_knots);
71 94 const int xnumpatch = unique_knots_x_.size() - 1;
72 94 const int ynumpatch = unique_knots_y_.size() - 1;
73 94 polynomial_degree_x_ = x_knots.size() - xyzw_cols;
74 94 polynomial_degree_y_ = y_knots.size() - xyzw_rows;
75 94 data_.resize(4 * (polynomial_degree_x_ * xnumpatch * polynomial_degree_y_ *
76 ynumpatch));
77 {
78 // Since its only for initialization, I do not care about speed.
79 // Here I look weather the data given is already in bezier form.
80
5/6
✓ Branch 1 taken 93 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 93 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 93 times.
✓ Branch 7 taken 1 times.
 94 if (unique_knots_x_.size() == 2 && unique_knots_y_.size() == 2) {
81
2/2
✓ Branch 0 taken 372 times.
✓ Branch 1 taken 93 times.
 465 for (int i = 0; i < 4; i++) {
82 372 Eigen::Matrix<double, -1, 1> tmp = Spl::Unroll(xyzw[i]);
83
3/4
✓ Branch 1 taken 5560 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 5560 times.
✓ Branch 6 taken 372 times.
 5932 for (int j = 0; j < tmp.rows(); j++) data_[j * 4 + i] = (tmp[j]);
84 372 }
85 } else {
86 // If not, I construct the dynamic projection (i.e. solve
87 // systems for
88 // the coeffs) and project to the superspace.
89
90 Eigen::SparseMatrix<double> phi = Spl::MakeProjection(
91 1 x_knots, y_knots, unique_knots_x_, unique_knots_y_,
92 1 polynomial_degree_x_, polynomial_degree_y_);
93
94
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (int i = 0; i < 4; i++) {
95 Eigen::Matrix<double, -1, 1> tmp =
96
4/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
 4 Spl::Unroll(xyzw[i]).transpose() * phi.transpose();
97
98
3/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 64 times.
✓ Branch 6 taken 4 times.
 68 for (int j = 0; j < tmp.rows(); j++) data_[j * 4 + i] = (tmp[j]);
99 4 }
100 1 }
101 }
102
103 94 return;
104 }
105 /**
106 * \brief Evaluate patch at given point.
107 *
108 * I look up the position in the knot vector, scale the input arguments,
109 * evaluate the 1D basis functions and sum over them with the control points
110 * from data.
111 *
112 * \param reference_point Point in reference domain.
113 * \return Point in physical domain.
114 */
115 80277 Eigen::Vector3d eval(const Eigen::Vector2d &reference_point) const {
116 const int x_location =
117
1/2
✓ Branch 1 taken 80277 times.
✗ Branch 2 not taken.
 80277 Spl::FindLocationInKnotVector(reference_point(0), unique_knots_x_);
118 const int y_location =
119
1/2
✓ Branch 1 taken 80277 times.
✗ Branch 2 not taken.
 80277 Spl::FindLocationInKnotVector(reference_point(1), unique_knots_y_);
120 80277 const int numy = (unique_knots_y_.size() - 1) * polynomial_degree_y_;
121 const double scaledx =
122
1/2
✓ Branch 2 taken 80277 times.
✗ Branch 3 not taken.
 80277 Spl::Rescale(reference_point(0), unique_knots_x_[x_location],
123 80277 unique_knots_x_[x_location + 1]);
124 const double scaledy =
125
1/2
✓ Branch 2 taken 80277 times.
✗ Branch 3 not taken.
 80277 Spl::Rescale(reference_point(1), unique_knots_y_[y_location],
126 80277 unique_knots_y_[y_location + 1]);
127
128
2/4
✓ Branch 0 taken 80277 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 80277 times.
✗ Branch 5 not taken.
 80277 double *xbasis = new double[polynomial_degree_x_];
129
2/4
✓ Branch 0 taken 80277 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 80277 times.
✗ Branch 5 not taken.
 80277 double *ybasis = new double[polynomial_degree_y_];
130
131
1/2
✓ Branch 1 taken 80277 times.
✗ Branch 2 not taken.
 80277 Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1,
132 xbasis, scaledx);
133
1/2
✓ Branch 1 taken 80277 times.
✗ Branch 2 not taken.
 80277 Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1,
134 ybasis, scaledy);
135
136 80277 double tmp[4] = {0., 0., 0., 0.};
137
138
2/2
✓ Branch 0 taken 167223 times.
✓ Branch 1 taken 80277 times.
 247500 for (int i = 0; i < polynomial_degree_x_; i++) {
139
2/2
✓ Branch 0 taken 368396 times.
✓ Branch 1 taken 167223 times.
 535619 for (int j = 0; j < polynomial_degree_y_; j++) {
140 368396 const double tpbasisval = xbasis[i] * ybasis[j];
141 368396 const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) +
142 368396 polynomial_degree_y_ * y_location + j);
143 #pragma omp simd
144
2/2
✓ Branch 1 taken 1473584 times.
✓ Branch 2 taken 368396 times.
 1841980 for (int k = 0; k < 4; k++) tmp[k] += data_[accs + k] * tpbasisval;
145 }
146 }
147
148
1/2
✓ Branch 0 taken 80277 times.
✗ Branch 1 not taken.
 80277 delete[] xbasis;
149
1/2
✓ Branch 0 taken 80277 times.
✗ Branch 1 not taken.
 80277 delete[] ybasis;
150
151
1/2
✓ Branch 1 taken 80277 times.
✗ Branch 2 not taken.
 80277 Eigen::Vector3d out(tmp[0], tmp[1], tmp[2]);
152 // Rescaling by the NRBS weight, i.e. projection to 3D from 4D hom
153
154
2/4
✓ Branch 1 taken 80277 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 80277 times.
✗ Branch 5 not taken.
 160554 return out / tmp[3];
155 }
156
157 /**
158 * \brief Evaluate Jacobian of the parametrization at given point.
159 *
160 * \param reference_point Point in reference domain.
161 * \return 3x2 Matrix with the Jacobian.
162 */
163 1313 Eigen::Matrix<double, 3, 2> evalJacobian(
164 const Eigen::Vector2d &reference_point) const {
165 const int x_location =
166
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Spl::FindLocationInKnotVector(reference_point(0), unique_knots_x_);
167 const int y_location =
168
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Spl::FindLocationInKnotVector(reference_point(1), unique_knots_y_);
169 1313 const int numy = (unique_knots_y_.size() - 1) * polynomial_degree_y_;
170 const double scaledx =
171
1/2
✓ Branch 2 taken 1313 times.
✗ Branch 3 not taken.
 1313 Spl::Rescale(reference_point(0), unique_knots_x_[x_location],
172 1313 unique_knots_x_[x_location + 1]);
173 const double scaledy =
174
1/2
✓ Branch 2 taken 1313 times.
✗ Branch 3 not taken.
 1313 Spl::Rescale(reference_point(1), unique_knots_y_[y_location],
175 1313 unique_knots_y_[y_location + 1]);
176
177
2/4
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1313 times.
✗ Branch 5 not taken.
 1313 double *xbasis = new double[polynomial_degree_x_];
178
2/4
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1313 times.
✗ Branch 5 not taken.
 1313 double *ybasis = new double[polynomial_degree_y_];
179
2/4
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1313 times.
✗ Branch 5 not taken.
 1313 double *xbasisD = new double[polynomial_degree_x_];
180
2/4
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1313 times.
✗ Branch 5 not taken.
 1313 double *ybasisD = new double[polynomial_degree_y_];
181
182
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1,
183 xbasis, scaledx);
184
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1,
185 ybasis, scaledy);
186
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_x_ - 1,
187 xbasisD, scaledx);
188
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_y_ - 1,
189 ybasisD, scaledy);
190
191 1313 double tmp[4] = {0., 0., 0., 0.};
192 1313 double tmpDx[4] = {0., 0., 0., 0.};
193 1313 double tmpDy[4] = {0., 0., 0., 0.};
194
195
2/2
✓ Branch 0 taken 3967 times.
✓ Branch 1 taken 1313 times.
 5280 for (int i = 0; i < polynomial_degree_x_; i++) {
196
2/2
✓ Branch 0 taken 15244 times.
✓ Branch 1 taken 3967 times.
 19211 for (int j = 0; j < polynomial_degree_y_; j++) {
197 15244 const double tpbasisval = xbasis[i] * ybasis[j];
198 15244 const double tpbasisvalDx = xbasisD[i] * ybasis[j];
199 15244 const double tpbasisvalDy = xbasis[i] * ybasisD[j];
200 15244 const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) +
201 15244 polynomial_degree_y_ * y_location + j);
202
203 // Here I add up the values of the basis functions in the dc
204 // basis
205 #pragma omp simd
206
2/2
✓ Branch 0 taken 60976 times.
✓ Branch 1 taken 15244 times.
 76220 for (int k = 0; k < 4; k++) {
207 60976 tmp[k] += data_[accs + k] * tpbasisval;
208 60976 tmpDx[k] += data_[accs + k] * tpbasisvalDx;
209 60976 tmpDy[k] += data_[accs + k] * tpbasisvalDy;
210 }
211 }
212 }
213
214
1/2
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
 1313 delete[] xbasis;
215
1/2
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
 1313 delete[] ybasis;
216
1/2
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
 1313 delete[] xbasisD;
217
1/2
✓ Branch 0 taken 1313 times.
✗ Branch 1 not taken.
 1313 delete[] ybasisD;
218
219
1/2
✓ Branch 1 taken 1313 times.
✗ Branch 2 not taken.
 1313 Eigen::Matrix<double, 3, 2> out;
220
221 // Eigen::Vector3d out;
222
223 1313 double bot = 1. / (tmp[3] * tmp[3]);
224
225 #pragma omp simd
226
2/2
✓ Branch 0 taken 3939 times.
✓ Branch 1 taken 1313 times.
 5252 for (int k = 0; k < 3; k++) {
227
1/2
✓ Branch 1 taken 3939 times.
✗ Branch 2 not taken.
 3939 out(k, 0) = (tmpDx[k] * tmp[3] - tmp[k] * tmpDx[3]) * bot;
228
1/2
✓ Branch 1 taken 3939 times.
✗ Branch 2 not taken.
 3939 out(k, 1) = (tmpDy[k] * tmp[3] - tmp[k] * tmpDy[3]) * bot;
229 }
230
231 2626 return out;
232 }
233
234 /**
235 * \brief Evaluate normal vector at given point.
236 *
237 * \param reference_point Point in reference domain.
238 * \return 3x1 normal vector.
239 */
240 1071 inline Eigen::Matrix<double, 3, 1> evalNormal(
241 const Eigen::Vector2d &reference_point) const {
242
1/2
✓ Branch 1 taken 1071 times.
✗ Branch 2 not taken.
 1071 Eigen::Matrix<double, 3, 2> jac = evalJacobian(reference_point);
243
3/6
✓ Branch 1 taken 1071 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1071 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1071 times.
✗ Branch 8 not taken.
 2142 return jac.col(0).cross(jac.col(1));
244 }
245
246 /**
247 * \brief Evaluate patch at given point.
248 *
249 * I look up the position in the knot vector, scale the input arguments,
250 * evaluate the 1D basis functions and sum over them with the control points
251 * from data.
252 *
253 * \param x Coordinate in reference domain.
254 * \param y Coordinate in reference domain.
255 * \return Point in physical domain.
256 */
257 77064 inline Eigen::Vector3d eval(double x, double y) const {
258
1/2
✓ Branch 2 taken 77064 times.
✗ Branch 3 not taken.
 154128 return eval(Eigen::Vector2d(x, y));
259 }
260 /**
261 * \brief Evaluate Jacobian of the parametrization at given point.
262 *
263 * \param x Coordinate in reference domain.
264 * \param y Coordinate in reference domain.
265 * \return 3x2 Matrix with the Jacobian.
266 */
267 inline Eigen::Matrix<double, 3, 2> evalJacobian(double x, double y) const {
268 return evalJacobian(Eigen::Vector2d(x, y));
269 }
270 /**
271 * \brief Evaluate normal vector at given point.
272 *
273 * \param x Coordinate in reference domain.
274 * \param y Coordinate in reference domain.
275 * \return 3x1 normal vector.
276 */
277 inline Eigen::Matrix<double, 3, 1> evalNormal(double x, double y) const {
278 return evalNormal(Eigen::Vector2d(x, y));
279 }
280
281 //
282 /**
283 * \brief Updates the surface point and returns the physical point and the
284 * derivatives there.
285 *
286 * This is a combination of eval und evalJacobian, to avoid duplication of
287 * work.
288 *
289 * \param srf_pt Pointer to the SurfacePoint which gets updated.
290 * \param ref_pt Point in reference domain with respect to the patch.
291 * \param w quadrature weight.
292 * \param xi Point in reference domain with respect to the element.
293 */
294 99824466 void updateSurfacePoint(SurfacePoint *srf_pt,
295 const Eigen::Vector2d &ref_pt, double w,
296 const Eigen::Vector2d &xi) const {
297 const int x_location =
298 99824466 Spl::FindLocationInKnotVector(ref_pt(0), unique_knots_x_);
299 const int y_location =
300 99824466 Spl::FindLocationInKnotVector(ref_pt(1), unique_knots_y_);
301 99824466 const int numy = (unique_knots_y_.size() - 1) * polynomial_degree_y_;
302 99824466 const double scaledx = Spl::Rescale(ref_pt(0), unique_knots_x_[x_location],
303 99824466 unique_knots_x_[x_location + 1]);
304 99824466 const double scaledy = Spl::Rescale(ref_pt(1), unique_knots_y_[y_location],
305 99824466 unique_knots_y_[y_location + 1]);
306
307 // TODO(Felix) Do not use variable-length arrays in accordance to Google
308 // Style guide
309 double *buffer =
310
1/2
✓ Branch 0 taken 99824466 times.
✗ Branch 1 not taken.
 99824466 new double[2 * (polynomial_degree_x_ + polynomial_degree_y_) + 12];
311
2/2
✓ Branch 0 taken 1197893592 times.
✓ Branch 1 taken 99824466 times.
 1297718058 for (int i = 0; i < 12; ++i) buffer[i] = 0;
312
313 99824466 double *tmp = buffer;
314 99824466 double *tmpDx = tmp + 4;
315 99824466 double *tmpDy = tmpDx + 4;
316 99824466 double *xbasis = tmpDy + 4;
317 99824466 double *ybasis = xbasis + polynomial_degree_x_;
318 99824466 double *xbasisD = ybasis + polynomial_degree_y_;
319 99824466 double *ybasisD = xbasisD + polynomial_degree_x_;
320
321 99824466 Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1,
322 xbasis, scaledx);
323 99824466 Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1,
324 ybasis, scaledy);
325 99824466 Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_x_ - 1,
326 xbasisD, scaledx);
327 99824466 Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_y_ - 1,
328 ybasisD, scaledy);
329
330
2/2
✓ Branch 0 taken 232534230 times.
✓ Branch 1 taken 99824466 times.
 332358696 for (int i = 0; i < polynomial_degree_x_; ++i) {
331
2/2
✓ Branch 0 taken 629494950 times.
✓ Branch 1 taken 232534230 times.
 862029180 for (int j = 0; j < polynomial_degree_y_; ++j) {
332 629494950 const double tpbasisval = xbasis[i] * ybasis[j];
333 629494950 const double tpbasisvalDx = xbasisD[i] * ybasis[j];
334 629494950 const double tpbasisvalDy = xbasis[i] * ybasisD[j];
335 629494950 const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) +
336 629494950 polynomial_degree_y_ * y_location + j);
337
338 // Here I add up the values of the basis functions in the dc
339 // basis
340
2/2
✓ Branch 0 taken 2517979800 times.
✓ Branch 1 taken 629494950 times.
 3147474750 for (int k = 0; k < 4; ++k) {
341 2517979800 tmp[k] += data_[accs + k] * tpbasisval;
342 2517979800 tmpDx[k] += data_[accs + k] * tpbasisvalDx;
343 2517979800 tmpDy[k] += data_[accs + k] * tpbasisvalDy;
344 }
345 }
346 }
347
348 99824466 const double bot = 1. / tmp[3];
349 99824466 const double botsqr = bot * bot;
350
351 99824466 (*srf_pt)(0) = xi(0);
352 99824466 (*srf_pt)(1) = xi(1);
353 99824466 (*srf_pt)(2) = w;
354 99824466 (*srf_pt)(3) = tmp[0] * bot;
355 99824466 (*srf_pt)(4) = tmp[1] * bot;
356 99824466 (*srf_pt)(5) = tmp[2] * bot;
357 99824466 (*srf_pt)(6) = (tmpDx[0] * tmp[3] - tmp[0] * tmpDx[3]) * botsqr;
358 99824466 (*srf_pt)(7) = (tmpDx[1] * tmp[3] - tmp[1] * tmpDx[3]) * botsqr;
359 99824466 (*srf_pt)(8) = (tmpDx[2] * tmp[3] - tmp[2] * tmpDx[3]) * botsqr;
360 99824466 (*srf_pt)(9) = (tmpDy[0] * tmp[3] - tmp[0] * tmpDy[3]) * botsqr;
361 99824466 (*srf_pt)(10) = (tmpDy[1] * tmp[3] - tmp[1] * tmpDy[3]) * botsqr;
362 99824466 (*srf_pt)(11) = (tmpDy[2] * tmp[3] - tmp[2] * tmpDy[3]) * botsqr;
363
1/2
✓ Branch 0 taken 99824466 times.
✗ Branch 1 not taken.
 99824466 delete[] buffer;
364 99824466 return;
365 }
366
367 //////////////////////////////////////////////////////////////////////////////
368 /// getter
369 //////////////////////////////////////////////////////////////////////////////
370
371 std::vector<double> data_; // Controllpoints in Bezier-Extracted Format.
372 int polynomial_degree_x_; // Degree in x
373 int polynomial_degree_y_; // Degree in y
374 std::vector<double>
375 unique_knots_x_; // The knot vectors, where each knot is unique
376 std::vector<double>
377 unique_knots_y_; // The knot vectors, where each knot is unique
378 };
379
380 /**
381 * \brief This function cuts a patch along internal knots, if any.
382 *
383 * \param patch Patch which gets checked and cut if there are internal knots.
384 * \return Vector of patches if it got cut, otherwise the vector is of size 1.
385 */
386 94 inline std::vector<Patch> PatchShredder(const Patch &patch) noexcept {
387 // Already a Bezier patch
388
5/6
✓ Branch 1 taken 93 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 93 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 93 times.
✓ Branch 7 taken 1 times.
 94 if (patch.unique_knots_y_.size() == 2 && patch.unique_knots_x_.size() == 2) {
389
2/2
✓ Branch 4 taken 93 times.
✓ Branch 5 taken 93 times.
 186 return {patch};
390 }
391
392 // number of subpatches in x and y directions
393 1 const int xchips = patch.unique_knots_x_.size() - 1;
394 1 const int ychips = patch.unique_knots_y_.size() - 1;
395
396 1 const int xp = patch.polynomial_degree_x_;
397 1 const int yp = patch.polynomial_degree_y_;
398 1 const int numy = ychips * yp;
399
400 1 std::vector<Patch> out(xchips * ychips);
401
402
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (int ix = 0; ix < xchips; ix++) {
403
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 for (int iy = 0; iy < ychips; iy++) {
404 4 const int index = ix * ychips + iy;
405
406 4 out[index].unique_knots_x_ = {0, 1};
407 4 out[index].unique_knots_y_ = {0, 1};
408 4 out[index].polynomial_degree_x_ = xp;
409 4 out[index].polynomial_degree_y_ = yp;
410 4 out[index].data_.reserve(xp * yp * 4);
411 }
412 }
413
414
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (int ix = 0; ix < xchips; ix++) {
415
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 for (int iy = 0; iy < ychips; iy++) {
416 4 const int index = ix * ychips + iy;
417
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 4 times.
 12 for (int jx = 0; jx < xp; jx++) {
418
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 8 times.
 24 for (int jy = 0; jy < yp; jy++) {
419 16 const int accs = 4 * (numy * (xp * ix + jx) + yp * iy + jy);
420
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 16 times.
 80 for (int k = 0; k < 4; k++) {
421 64 out[index].data_.push_back(patch.data_[accs + k]);
422 }
423 }
424 }
425 }
426 }
427
428 1 return out;
429 1 }
430
431 /**
432 * \brief This function cuts all patches along internal knots, if any.
433 *
434 * \param patches Vector with patches.
435 * \return Vector of patches if it got cut.
436 */
437 27 inline std::vector<Patch> PatchShredder(
438 const std::vector<Patch> &patches) noexcept {
439 27 std::vector<Patch> out;
440 27 const int input_size = patches.size();
441
442
2/2
✓ Branch 0 taken 94 times.
✓ Branch 1 taken 27 times.
 121 for (int i = 0; i < input_size; i++) {
443 94 std::vector<Patch> tmp = PatchShredder(patches[i]);
444
445 94 const int tmp_size = tmp.size();
446
447
2/2
✓ Branch 2 taken 97 times.
✓ Branch 3 taken 94 times.
 191 for (int j = 0; j < tmp_size; j++) out.push_back(tmp[j]);
448 94 }
449
450 27 out.shrink_to_fit();
451
452 27 return out;
453 }
454
455 } // namespace Bembel
456
457 #endif // BEMBEL_SRC_GEOMETRY_PATCH_HPP_
458