GCC Code Coverage Report

Directory: Bembel/src/
File: ClusterTree/ElementTree.hpp
Date: 2024-12-18 07:36:36
Exec Total Coverage
Lines: 224 234 95.7%
Functions: 19 19 100.0%
Branches: 182 290 62.8%
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_CLUSTERTREE_ELEMENTTREE_HPP_
13 #define BEMBEL_SRC_CLUSTERTREE_ELEMENTTREE_HPP_
14
15 namespace Bembel {
16 /**
17 * \ingroup ClusterTree
18 * \brief This class organizes an element structure on a Geometry object and
19 * handles refinement.
20 *
21 * This class implements a tree which realizes the refinement of the geometry.
22 * On the first level of refinement the leafs/elements are all patches of the
23 * geometry. Currently, the construction of the tree allows uniform refinement.
24 * In the refinement process each element get four sons and the neighborhood
25 * relations get updated.
26 *
27 * The heart of the ElementTree is the leaf iterator which implements a Morton
28 * Z-curve. If an element get refined it gets replaced by it four sons element.
29 * This core routine can handle adaptive refinement but with some limitations in
30 * the resolution of neighborhood relations.
31 */
32 class ElementTree {
33 public:
34 //////////////////////////////////////////////////////////////////////////////
35 /**
36 * \brief Deleted copy constructor for the ElementTree class.
37 *
38 * This copy constructor is explicitly deleted to prevent copying of
39 * ElementTree objects.
40 */
41 ElementTree(const ElementTree &) = delete;
42
43 /**
44 * \brief Deleted move constructor for the ElementTree class.
45 *
46 * This move constructor is explicitly deleted to prevent moving of
47 * ElementTree objects.
48 */
49 ElementTree(ElementTree &&) = delete;
50
51 /**
52 * \brief Deleted copy assignment operator for the ElementTree class.
53 *
54 * This copy assignment operator is explicitly deleted to prevent copying of
55 * ElementTree objects.
56 *
57 * \return A reference to the updated ElementTree object.
58 */
59 ElementTree &operator=(const ElementTree &) = delete;
60
61 /**
62 * \brief Deleted move assignment operator for the ElementTree class.
63 *
64 * This move assignment operator is explicitly deleted to prevent moving of
65 * ElementTree objects.
66 *
67 * \return A reference to the updated ElementTree object.
68 */
69 ElementTree &operator=(ElementTree &&) = delete;
70 //////////////////////////////////////////////////////////////////////////////
71 /// constructors
72 //////////////////////////////////////////////////////////////////////////////
73 /**
74 * \brief Default constructor for the ElementTree class.
75 */
76 168 ElementTree() {}
77 /**
78 * \brief Explicit constructor for the ElementTree class.
79 *
80 * This constructor initializes an ElementTree object for the provided
81 * geometry and maximum level. Currently this constructor is implemented for
82 * uniform refinement.
83 *
84 * \param g The geometry object defining the patches.
85 * \param max_level (optional) The maximum level of the ElementTree. Default
86 * is 0.
87 */
88 explicit ElementTree(const Geometry &g, unsigned int max_level = 0) {
89 init_ElementTree(g, max_level);
90 }
91 //////////////////////////////////////////////////////////////////////////////
92 /// init
93 //////////////////////////////////////////////////////////////////////////////
94 /**
95 * \brief This function initializes the ElementTree
96 *
97 * First Each Patch becomes an ElementTreeNode. After resolving the patch
98 * neighborhood relations each patch gets refined to the common maximum level
99 * of refinement.
100 *
101 * \param g The geometry object defining the patches.
102 * \param max_level (optional) The maximum level of the ElementTree. Default
103 * is 0.
104 */
105 168 void init_ElementTree(const Geometry &g, unsigned int max_level) {
106 // initialise the data fields
107 168 geometry_ = g.get_geometry_ptr();
108 168 max_level_ = max_level;
109 168 number_of_patches_ = geometry_->size();
110 // create the patches and set up the topology
111 {
112 168 std::vector<Eigen::Vector3d> uniquePts;
113 168 std::vector<ElementTreeNode *> patches;
114
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 168 Eigen::Vector3d v;
115 168 number_of_points_ = 0;
116 168 number_of_elements_ = number_of_patches_;
117 168 ElementTreeNode &root = root_;
118
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 168 root.sons_.resize(number_of_patches_);
119
2/2
✓ Branch 0 taken 729 times.
✓ Branch 1 taken 168 times.
 897 for (auto i = 0; i < number_of_patches_; ++i) {
120
1/2
✓ Branch 2 taken 729 times.
✗ Branch 3 not taken.
 729 root.sons_[i].adjcents_ = std::vector<ElementTreeNode *>(4, nullptr);
121
1/2
✓ Branch 2 taken 729 times.
✗ Branch 3 not taken.
 729 root.sons_[i].vertices_.resize(4);
122 729 root.sons_[i].id_ = i;
123 729 root.sons_[i].level_ = 0;
124 729 root.sons_[i].patch_ = i;
125 // add linked list structure to the panels
126
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 561 times.
 729 if (i == 0) {
127 168 root.sons_[0].prev_ = nullptr;
128 168 pfirst_ = std::addressof(root.sons_[0]);
129 } else {
130 561 root.sons_[i].prev_ = std::addressof(root.sons_[i - 1]);
131 }
132
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 561 times.
 729 if (i == number_of_patches_ - 1) {
133 168 root.sons_[i].next_ = nullptr;
134 168 plast_ = std::addressof(root.sons_[i]);
135 } else {
136 561 root.sons_[i].next_ = std::addressof(root.sons_[i + 1]);
137 }
138 // get images of the four corners of the unit square under the
139 // diffeomorphism geo[i]
140
2/2
✓ Branch 0 taken 2916 times.
✓ Branch 1 taken 729 times.
 3645 for (auto j = 0; j < 4; ++j) {
141 2916 v = (*geometry_)[i].eval(Constants::corners[0][j],
142
1/2
✓ Branch 1 taken 2916 times.
✗ Branch 2 not taken.
 2916 Constants::corners[1][j]);
143 2916 unsigned int index = 0;
144
2/2
✓ Branch 1 taken 12174 times.
✓ Branch 2 taken 1740 times.
 13914 for (; index < uniquePts.size(); ++index)
145
4/6
✓ Branch 2 taken 12174 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12174 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1176 times.
✓ Branch 8 taken 10998 times.
 12174 if ((uniquePts[index] - v).norm() < Constants::pt_comp_tolerance)
146 1176 break;
147
2/2
✓ Branch 1 taken 1176 times.
✓ Branch 2 taken 1740 times.
 2916 if (index != uniquePts.size()) {
148 1176 root.sons_[i].vertices_[j] = index;
149 } else {
150
1/2
✓ Branch 1 taken 1740 times.
✗ Branch 2 not taken.
 1740 uniquePts.push_back(v);
151 1740 root.sons_[i].vertices_[j] = number_of_points_;
152 1740 ++number_of_points_;
153 }
154 }
155
1/2
✓ Branch 3 taken 729 times.
✗ Branch 4 not taken.
 729 patches.push_back(std::addressof(root.sons_[i]));
156 }
157
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 168 updateTopology(patches);
158 168 }
159
2/2
✓ Branch 1 taken 241 times.
✓ Branch 2 taken 168 times.
 409 for (auto i = 0; i < max_level; ++i) refineUniformly();
160 168 return;
161 }
162 //////////////////////////////////////////////////////////////////////////////
163 /**
164 * \brief Refine the given element or it's sons.
165 *
166 * This function refines an element recursively. If the given element is
167 * already refined than iterate all sons and refine them and so on.
168 */
169 7801 void refineUniformly_recursion(ElementTreeNode &el) {
170
2/2
✓ Branch 1 taken 1865 times.
✓ Branch 2 taken 5936 times.
 7801 if (el.sons_.size()) {
171
2/2
✓ Branch 1 taken 7560 times.
✓ Branch 2 taken 1865 times.
 9425 for (auto i = 0; i < el.sons_.size(); ++i)
172 7560 refineUniformly_recursion(el.sons_[i]);
173 } else {
174 5936 refineLeaf(el);
175 }
176 7801 return;
177 }
178 //////////////////////////////////////////////////////////////////////////////
179 /**
180 * \brief Refine all patches uniformly
181 */
182 241 void refineUniformly() {
183 241 refineUniformly_recursion(root_);
184 241 return;
185 }
186 //////////////////////////////////////////////////////////////////////////////
187 /**
188 * \brief Refine a given patch.
189 */
190 void refinePatch(int patch) {
191 refineUniformly_recursion(root_.sons_[patch]);
192 return;
193 }
194 //////////////////////////////////////////////////////////////////////////////
195 /**
196 * \brief Return the coordinates of all points of the elements.
197 *
198 * This function iterates all elements and returns the coordinates of the
199 * vertices.
200 *
201 * \param idct Pointer to an Eigen::VectorXi to count how often a vertex is
202 * part of all elements.
203 *
204 * \return A 3xN Matrix where N is the number of all vertices in the geometry.
205 */
206 168 Eigen::MatrixXd generatePointList(Eigen::VectorXi *idct = nullptr) const {
207
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 168 Eigen::MatrixXd pts(3, number_of_points_);
208
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 168 times.
 168 if (idct != nullptr) {
209 ✗ idct->resize(number_of_points_);
210 ✗ idct->setZero();
211 }
212
2/2
✓ Branch 4 taken 18537 times.
✓ Branch 5 taken 168 times.
 18705 for (auto it = pbegin(); it != pend(); ++it) {
213 18537 double h = double(1) / double(1 << it->level_);
214
1/2
✓ Branch 3 taken 18537 times.
✗ Branch 4 not taken.
 18537 pts.col(it->vertices_[0]) =
215
4/8
✓ Branch 5 taken 18537 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18537 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18537 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18537 times.
✗ Branch 16 not taken.
 37074 (*geometry_)[it->patch_].eval(it->llc_(0), it->llc_(1));
216
1/2
✓ Branch 3 taken 18537 times.
✗ Branch 4 not taken.
 18537 pts.col(it->vertices_[1]) =
217
4/8
✓ Branch 5 taken 18537 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18537 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18537 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18537 times.
✗ Branch 16 not taken.
 37074 (*geometry_)[it->patch_].eval(it->llc_(0) + h, it->llc_(1));
218
1/2
✓ Branch 3 taken 18537 times.
✗ Branch 4 not taken.
 18537 pts.col(it->vertices_[2]) =
219
4/8
✓ Branch 5 taken 18537 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18537 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18537 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18537 times.
✗ Branch 16 not taken.
 37074 (*geometry_)[it->patch_].eval(it->llc_(0) + h, it->llc_(1) + h);
220
1/2
✓ Branch 3 taken 18537 times.
✗ Branch 4 not taken.
 18537 pts.col(it->vertices_[3]) =
221
4/8
✓ Branch 5 taken 18537 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18537 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18537 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18537 times.
✗ Branch 16 not taken.
 37074 (*geometry_)[it->patch_].eval(it->llc_(0), it->llc_(1) + h);
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18537 times.
 18537 if (idct != nullptr) {
223 ✗ ++((*idct)(it->vertices_[0]));
224 ✗ ++((*idct)(it->vertices_[1]));
225 ✗ ++((*idct)(it->vertices_[2]));
226 ✗ ++((*idct)(it->vertices_[3]));
227 }
228 }
229 168 return pts;
230 ✗ }
231 //////////////////////////////////////////////////////////////////////////////
232 /**
233 * \brief Return list of elements containing the indices of the vertices.
234 *
235 * This function stores all indices of an element in the column of the
236 * returned matrix.
237 *
238 * \return A 4xN integer Matrix where N is the number elements.
239 */
240 Eigen::MatrixXi generateElementList() const {
241 unsigned int elID = 0;
242 Eigen::MatrixXi retval(4, number_of_elements_);
243 for (auto it = pbegin(); it != pend(); ++it) {
244 retval.col(elID) << it->vertices_[0], it->vertices_[1], it->vertices_[2],
245 it->vertices_[3];
246 ++elID;
247 }
248 return retval;
249 }
250 //////////////////////////////////////////////////////////////////////////////
251 /// getter
252 //////////////////////////////////////////////////////////////////////////////
253 /**
254 * \brief Return number of points in the ElementTree.
255 *
256 * \return Number of points.
257 */
258 int get_number_of_points() const { return number_of_points_; }
259 /**
260 * \brief Return number of elements in the ElementTree.
261 *
262 * \return Number of elements.
263 */
264 40 int get_number_of_elements() const { return number_of_elements_; }
265 /**
266 * \brief Return maximum level of refinement.
267 *
268 * \return Level of refinement.
269 */
270 2567942 int get_max_level() const { return max_level_; }
271 /**
272 * \brief Return reference to the root ElementTreeNode.
273 *
274 * \return Reference to the root ElementTreeNode.
275 */
276 ElementTreeNode &root() { return root_; }
277 /**
278 * \brief Return const reference to the root ElementTreeNode.
279 *
280 * \return Const Reference to the root ElementTreeNode.
281 */
282 15 const ElementTreeNode &root() const { return root_; }
283 /**
284 * \brief Return const reference to the Geometry.
285 *
286 * \return Const Reference to the PatchVector.
287 */
288 99825513 const PatchVector &get_geometry() const { return *geometry_; }
289 //////////////////////////////////////////////////////////////////////////////
290 /// iterators
291 //////////////////////////////////////////////////////////////////////////////
292 /**
293 * \brief Returns an iterator to the beginning of the sequence represented by
294 * the leafs as ElementTreeNodes of the ElementTree.
295 *
296 * \return Returns a ElementTreeNode::const_iterator object.
297 */
298 4585 ElementTreeNode::const_iterator pbegin() const {
299 4585 return ElementTreeNode::const_iterator(pfirst_);
300 }
301 /**
302 * \brief Returns an iterator one past the end of the sequence represented by
303 * the leafs as ElementTreeNodes of the ElementTree.
304 *
305 * \return Returns a ElementTreeNode::const_iterator object.
306 */
307 69734 ElementTreeNode::const_iterator pend() const {
308 69734 return ElementTreeNode::const_iterator(plast_->next_);
309 }
310 /**
311 * \brief Returns an iterator to the beginning of the sequence represented by
312 * the leafs as ElementTreeNodes of the ElementTree.
313 *
314 * \return Returns a ElementTreeNode::const_iterator object.
315 */
316 4271 ElementTreeNode::const_iterator cpbegin() const { return pbegin(); }
317 /**
318 * \brief Returns an iterator one past the end of the sequence represented by
319 * the leafs as ElementTreeNodes of the ElementTree.
320 *
321 * \return Returns a ElementTreeNode::const_iterator object.
322 */
323 36863 ElementTreeNode::const_iterator cpend() const { return pend(); }
324 /**
325 * \brief Returns a cluster iterator to the beginning of the sequence
326 * represented by the the given ElementTreeNode.
327 *
328 * \param cl Const reference to the ElementTreeNode to start the iterator.
329 * \return Returns a ElementTreeNode::const_iterator object.
330 */
331 ElementTreeNode::const_iterator cluster_begin(
332 const ElementTreeNode &cl) const {
333 return cl.cbegin();
334 }
335 /**
336 * \brief Returns a cluster iterator one past the end of the sequence
337 * represented by the the given ElementTreeNode.
338 *
339 * \param cl Const reference to the ElementTreeNode to start the iterator.
340 * \return Returns a ElementTreeNode::const_iterator object.
341 */
342 ElementTreeNode::const_iterator cluster_end(const ElementTreeNode &cl) const {
343 return cl.cend();
344 }
345 //////////////////////////////////////////////////////////////////////////////
346 /**
347 * \brief Computes enclosing balls surrounding all elements.
348 *
349 * This functions sets the parameters midpoint_ radius_ of the ElementTreeNode
350 * objects stored in the tree.
351 *
352 * \return The point list as 3xN matrix with N the number of points.
353 */
354 168 Eigen::MatrixXd computeElementEnclosings() {
355 // compute point list
356 168 Eigen::MatrixXd P = generatePointList();
357
2/2
✓ Branch 1 taken 729 times.
✓ Branch 2 taken 168 times.
 897 for (auto i = 0; i < root_.sons_.size(); ++i)
358
1/2
✓ Branch 2 taken 729 times.
✗ Branch 3 not taken.
 729 computeElementEnclosings_recursion(root_.sons_[i], P);
359 168 return P;
360 ✗ }
361 //////////////////////////////////////////////////////////////////////////////
362 /**
363 * \brief Computes enclosing balls of one element and all its sons.
364 *
365 * \param el ElementTreeNode to start the recursion.
366 * \param P Point list of the vertices.
367 */
368 24473 void computeElementEnclosings_recursion(ElementTreeNode &el,
369 const Eigen::MatrixXd &P) {
370
2/4
✓ Branch 1 taken 24473 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24473 times.
✗ Branch 5 not taken.
 24473 Eigen::Vector3d mp1, mp2;
371 double r1, r2;
372 // compute point list
373
2/2
✓ Branch 1 taken 18537 times.
✓ Branch 2 taken 5936 times.
 24473 if (!el.sons_.size()) {
374 // assign enclosing balls to leafs
375
4/8
✓ Branch 1 taken 18537 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 18537 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 18537 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 18537 times.
✗ Branch 12 not taken.
 18537 util::computeEnclosingBall(&mp1, &r1, P.col(el.vertices_[0]), 0,
376
1/2
✓ Branch 2 taken 18537 times.
✗ Branch 3 not taken.
 18537 P.col(el.vertices_[2]), 0);
377
4/8
✓ Branch 1 taken 18537 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 18537 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 18537 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 18537 times.
✗ Branch 12 not taken.
 18537 util::computeEnclosingBall(&mp2, &r2, P.col(el.vertices_[1]), 0,
378
1/2
✓ Branch 2 taken 18537 times.
✗ Branch 3 not taken.
 18537 P.col(el.vertices_[3]), 0);
379
1/2
✓ Branch 1 taken 18537 times.
✗ Branch 2 not taken.
 18537 util::computeEnclosingBall(&(el.midpoint_), &(el.radius_), mp1, r1, mp2,
380 r2);
381 } else {
382 // handle the four(!!!) children
383
2/2
✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
 29680 for (auto i = 0; i < 4; ++i)
384
1/2
✓ Branch 2 taken 23744 times.
✗ Branch 3 not taken.
 23744 computeElementEnclosings_recursion(el.sons_[i], P);
385 // assign enclosing balls to fathers bottom up
386
1/2
✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
 5936 util::computeEnclosingBall(&mp1, &r1, el.sons_[0].midpoint_,
387 5936 el.sons_[0].radius_, el.sons_[2].midpoint_,
388 5936 el.sons_[2].radius_);
389
1/2
✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
 5936 util::computeEnclosingBall(&mp2, &r2, el.sons_[1].midpoint_,
390 5936 el.sons_[1].radius_, el.sons_[3].midpoint_,
391 5936 el.sons_[3].radius_);
392
1/2
✓ Branch 1 taken 5936 times.
✗ Branch 2 not taken.
 5936 util::computeEnclosingBall(&(el.midpoint_), &(el.radius_), mp1, r1, mp2,
393 r2);
394 }
395 48946 return;
396 }
397 //////////////////////////////////////////////////////////////////////////////
398 /**
399 * \brief Prints all Elements of the Tree.
400 */
401 void printPanels() const {
402 auto i = 0;
403 for (auto it = pbegin(); it != pend(); ++it) {
404 std::cout << "element[" << i++ << "] = ";
405 it->print();
406 }
407 return;
408 }
409 //////////////////////////////////////////////////////////////////////////////
410 /// other Stuff
411 //////////////////////////////////////////////////////////////////////////////
412 /**
413 * \brief Return a matrix with all midpoints of the elements.
414 *
415 * \return The midpoint list as 3xN matrix with N the number of points.
416 */
417 Eigen::MatrixXd generateMidpointList() const {
418 Eigen::MatrixXd retval(3, number_of_elements_);
419 unsigned int i = 0;
420 for (auto it = pbegin(); it != pend(); ++it) {
421 retval.col(i) = it->midpoint_;
422 ++i;
423 }
424 return retval;
425 }
426 //////////////////////////////////////////////////////////////////////////////
427 /**
428 * \brief Return a matrix with all radii of the element enclosing.
429 *
430 * \return A vector containing the radii of the element enclosing.
431 */
432 Eigen::MatrixXd generateRadiusList() const {
433 Eigen::VectorXd retval(number_of_elements_);
434 unsigned int i = 0;
435 for (auto it = pbegin(); it != pend(); ++it) {
436 retval(i) = it->radius_;
437 ++i;
438 }
439 return retval;
440 }
441 //////////////////////////////////////////////////////////////////////////////
442 /**
443 * \brief Return a vector with computed global indices of the elements.
444 *
445 * \return A vector with global indices of the elements.
446 */
447 Eigen::VectorXi generateElementLabels() const {
448 Eigen::VectorXi retval(number_of_elements_);
449 unsigned int i = 0;
450 for (auto it = pbegin(); it != pend(); ++it) {
451 retval(i) = compute_global_id(*it);
452 ++i;
453 }
454 return retval;
455 }
456 //////////////////////////////////////////////////////////////////////////////
457 /**
458 * \brief Generate list of labels if elements are on the patch boundary or at
459 * the boundary of the geometry.
460 *
461 * The value at the index of the element is 1 if the element is at the patch
462 * boundary and there is a neighbor patch. If there is no neighbor patch then
463 * the value is -1.
464 *
465 * \return A vector of integers of size N with N the number of elements.
466 */
467 Eigen::VectorXi generatePatchBoundaryLabels() const {
468 Eigen::VectorXi retval(number_of_elements_);
469 retval.setZero();
470 unsigned int i = 0;
471 for (auto it = pbegin(); it != pend(); ++it) {
472 for (auto j = 0; j < 4; ++j)
473 if (it->adjcents_[j] == nullptr) {
474 retval[i] = -1;
475 break;
476 } else if (it->adjcents_[j]->patch_ != it->patch_) {
477 ++(retval[i]);
478 }
479 ++i;
480 }
481 return retval;
482 }
483 //////////////////////////////////////////////////////////////////////////////
484 /**
485 * \brief Generate list of labels if elements contained within a given patch.
486 *
487 * The value at the index of the element is 1 if the element is contained
488 * within the given patch. Otherwise its zero.
489 *
490 * \return A vector of integers of size N with N the number of elements.
491 */
492 Eigen::VectorXi identifyPatch(unsigned int pn) const {
493 Eigen::VectorXi retval(number_of_elements_);
494 retval.setZero();
495 assert(pn < number_of_patches_);
496 unsigned int i = 0;
497 for (auto it = pbegin(); it != pend(); ++it) {
498 retval[i] = it->patch_ == pn ? 1 : 0;
499 ++i;
500 }
501 return retval;
502 }
503 /**
504 * \brief Resolves neighborhood relations of the patches.
505 *
506 * \return A vector where each entry defines a patch interface or boundary.
507 * The entries correspond to [patchIndex1, edgeCase1, patchIndex2, edgeCase2].
508 */
509 //////////////////////////////////////////////////////////////////////////////
510 312 std::vector<std::array<int, 4>> patchTopologyInfo() const {
511 312 std::vector<std::array<int, 4>> retval;
512
2/2
✓ Branch 4 taken 1434 times.
✓ Branch 5 taken 312 times.
 1746 for (auto it = root_.sons_.begin(); it != root_.sons_.end(); ++it) {
513
2/2
✓ Branch 0 taken 5736 times.
✓ Branch 1 taken 1434 times.
 7170 for (auto j = 0; j < 4; ++j) {
514 // do we have a neighbour?
515
2/2
✓ Branch 2 taken 2496 times.
✓ Branch 3 taken 3240 times.
 5736 if (it->adjcents_[j] != nullptr) {
516 2496 const ElementTreeNode &cur_neighbour = *(it->adjcents_[j]);
517 // add the edge only if it->id_ < neighbour->id_
518 // otherwise the edge has already been added
519
2/2
✓ Branch 1 taken 1248 times.
✓ Branch 2 taken 1248 times.
 2496 if (it->id_ < cur_neighbour.id_) {
520 1248 int k = 0;
521
1/2
✓ Branch 0 taken 3106 times.
✗ Branch 1 not taken.
 3106 for (; k < 4; ++k)
522
4/4
✓ Branch 1 taken 1558 times.
✓ Branch 2 taken 1548 times.
✓ Branch 3 taken 1248 times.
✓ Branch 4 taken 1858 times.
 4664 if (cur_neighbour.adjcents_[k] != nullptr &&
523
2/2
✓ Branch 2 taken 1248 times.
✓ Branch 3 taken 310 times.
 1558 cur_neighbour.adjcents_[k]->id_ == it->id_)
524 1248 break;
525
1/2
✓ Branch 2 taken 1248 times.
✗ Branch 3 not taken.
 1248 retval.push_back({it->id_, cur_neighbour.id_, j, k});
526 }
527 } else {
528
1/2
✓ Branch 2 taken 3240 times.
✗ Branch 3 not taken.
 3240 retval.push_back({it->id_, -1, j, -1});
529 }
530 }
531 }
532 312 return retval;
533 ✗ }
534 /**
535 * \brief The ordering of elements in the element tree does not correspond to
536 * the element order underlying the coefficient vector. This reordering can be
537 * computed for look ups by this function.
538 *
539 * Limitation to the uniform case!
540 *
541 * \return Vector with the tensor product index of the elements.
542 */
543 //////////////////////////////////////////////////////////////////////////////
544 138 std::vector<int> computeReorderingVector() const {
545
1/2
✓ Branch 2 taken 138 times.
✗ Branch 3 not taken.
 138 std::vector<int> out(number_of_elements_);
546
2/2
✓ Branch 3 taken 14028 times.
✓ Branch 4 taken 138 times.
 14166 for (auto it = pbegin(); it != pend(); ++it) {
547 14028 const double h = it->get_h();
548 const Eigen::Vector2d ref_mid =
549
3/6
✓ Branch 1 taken 14028 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 14028 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 14028 times.
✗ Branch 9 not taken.
 14028 it->llc_ + Eigen::Vector2d(.5 * h, .5 * h);
550
1/2
✓ Branch 1 taken 14028 times.
✗ Branch 2 not taken.
 14028 const int x_idx = std::floor(ref_mid(0) / h);
551
1/2
✓ Branch 1 taken 14028 times.
✗ Branch 2 not taken.
 14028 const int y_idx = std::floor(ref_mid(1) / h);
552 14028 const int num_in_one_dir = 1 << max_level_;
553 14028 const int num_per_patch = num_in_one_dir * num_in_one_dir;
554 const int tp_idx =
555 14028 it->patch_ * num_per_patch + y_idx * num_in_one_dir + x_idx;
556 14028 out[tp_idx] = it->id_;
557 }
558 138 return out;
559 ✗ }
560 //////////////////////////////////////////////////////////////////////////////
561 size_t compute_global_id(const ElementTreeNode &el) const {
562 return number_of_patches_ *
563 (((1 << el.level_) * (1 << el.level_) - 1) / 3) +
564 el.id_;
565 }
566 //////////////////////////////////////////////////////////////////////////////
567 /// private members
568 //////////////////////////////////////////////////////////////////////////////
569 private:
570 //////////////////////////////////////////////////////////////////////////////
571 /// methods
572 //////////////////////////////////////////////////////////////////////////////
573 template <typename T>
574 struct isEqual {
575 bool operator()(const T &v1, const T &v2) const {
576 return ((v1 - v2).norm() < Constants::pt_comp_tolerance);
577 }
578 };
579 /**
580 * \brief function to set up the local topology, i.e. the adjacents
581 * of a refined element
582 */
583 //////////////////////////////////////////////////////////////////////////////
584 6104 void updateTopology(const std::vector<ElementTreeNode *> &elements) {
585 6104 std::map<std::array<int, 2>, ElementTreeNode *> edges;
586 std::array<int, 2> e1, e2;
587 // generate edge list for all elements in question
588
2/2
✓ Branch 1 taken 43545 times.
✓ Branch 2 taken 6104 times.
 49649 for (auto i = 0; i < elements.size(); ++i)
589
2/2
✓ Branch 0 taken 174180 times.
✓ Branch 1 taken 43545 times.
 217725 for (auto j = 0; j < 4; ++j) {
590 // compute a unique id for each edge
591 174180 const int v1 = elements[i]->vertices_[j];
592 174180 const int v2 = elements[i]->vertices_[(j + 1) % 4];
593
4/4
✓ Branch 0 taken 87249 times.
✓ Branch 1 taken 86931 times.
✓ Branch 2 taken 87249 times.
✓ Branch 3 taken 86931 times.
 174180 e1 = {v1 < v2 ? v1 : v2, v1 < v2 ? v2 : v1};
594 // perferm a look up if the edge is already existing.
595
1/2
✓ Branch 1 taken 174180 times.
✗ Branch 2 not taken.
 174180 auto it = edges.find(e1);
596 // if so, identify the two neighbours
597
2/2
✓ Branch 2 taken 53064 times.
✓ Branch 3 taken 121116 times.
 174180 if (it != edges.end()) {
598 53064 elements[i]->adjcents_[j] = it->second;
599 // now find the edge also in the patch that added it to the set
600
1/2
✓ Branch 0 taken 156511 times.
✗ Branch 1 not taken.
 156511 for (auto k = 0; k < 4; ++k) {
601 156511 const int v3 = it->second->vertices_[k];
602 156511 const int v4 = it->second->vertices_[(k + 1) % 4];
603
4/4
✓ Branch 0 taken 86679 times.
✓ Branch 1 taken 69832 times.
✓ Branch 2 taken 86679 times.
✓ Branch 3 taken 69832 times.
 156511 e2 = {v3 < v4 ? v3 : v4, v3 < v4 ? v4 : v3};
604
3/4
✓ Branch 1 taken 156511 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 53064 times.
✓ Branch 4 taken 103447 times.
 156511 if (e1 == e2) {
605 53064 it->second->adjcents_[k] = elements[i];
606 53064 break;
607 }
608 }
609 // otherwise add the edge to the list
610 } else {
611
2/4
✓ Branch 2 taken 121116 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 121116 times.
✗ Branch 6 not taken.
 121116 edges.insert(std::make_pair(e1, elements[i]));
612 }
613 }
614 12208 return;
615 6104 }
616 //////////////////////////////////////////////////////////////////////////////
617 /**
618 * \brief This function refines the given ElementTreeNode.
619 *
620 * The function introduces 4 new elements and takes care of the newly
621 * introduces vertices. Furthermore, all the neighborhood relation of all
622 * surrounding elements are resolved.
623 */
624 5936 void refineLeaf(ElementTreeNode &cur_el) {
625 // check if we have actually a panel
626
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5936 times.
 5936 if (cur_el.sons_.size()) return;
627
1/2
✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
 5936 std::vector<int> ptIds(5);
628
1/2
✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
 5936 std::vector<int> refNeighbours(4);
629 5936 std::vector<ElementTreeNode *> elements;
630 // determine position of p with respect to its neighbours
631
2/2
✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
 29680 for (auto i = 0; i < 4; ++i) {
632 // is there a neighbour?
633
2/2
✓ Branch 1 taken 19072 times.
✓ Branch 2 taken 4672 times.
 23744 if (cur_el.adjcents_[i] != nullptr) {
634
1/2
✓ Branch 0 taken 47680 times.
✗ Branch 1 not taken.
 47680 for (auto j = 0; j < 4; ++j)
635
2/2
✓ Branch 3 taken 19072 times.
✓ Branch 4 taken 28608 times.
 47680 if (cur_el.adjcents_[i]->adjcents_[j] == std::addressof(cur_el)) {
636 19072 refNeighbours[i] = j;
637 19072 break;
638 }
639 } else {
640 4672 refNeighbours[i] = -1;
641 }
642 }
643 // determine new points
644
2/2
✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
 29680 for (auto i = 0; i < 4; ++i) {
645 // is there a neighbour?
646
2/2
✓ Branch 1 taken 19072 times.
✓ Branch 2 taken 4672 times.
 23744 if (cur_el.adjcents_[i] != nullptr) {
647 19072 ElementTreeNode &ref_cur_neighbour = *(cur_el.adjcents_[i]);
648 // is the neighbour already refined?
649
2/2
✓ Branch 1 taken 9536 times.
✓ Branch 2 taken 9536 times.
 19072 if (ref_cur_neighbour.sons_.size()) {
650 // this is the midpoint of the shared edge
651 19072 ptIds[i] = ref_cur_neighbour.sons_[refNeighbours[i]]
652 19072 .vertices_[(refNeighbours[i] + 1) % 4];
653 // these are the two elements adjacent to the edge
654 9536 elements.push_back(
655
1/2
✓ Branch 4 taken 9536 times.
✗ Branch 5 not taken.
 9536 std::addressof(ref_cur_neighbour.sons_[refNeighbours[i]]));
656
1/2
✓ Branch 2 taken 9536 times.
✗ Branch 3 not taken.
 9536 elements.push_back(std::addressof(
657 9536 ref_cur_neighbour.sons_[(refNeighbours[i] + 1) % 4]));
658 } else {
659 // otherwise add the point id to the tree
660 9536 ptIds[i] = number_of_points_;
661 9536 ++number_of_points_;
662 }
663 } else {
664 // otherwise add the point id to the tree
665 4672 ptIds[i] = number_of_points_;
666 4672 ++number_of_points_;
667 }
668 }
669 // add midpoint of the current element, which is always a new point
670 5936 ptIds[4] = number_of_points_;
671 5936 ++number_of_points_;
672 // set up new elements
673
1/2
✓ Branch 1 taken 5936 times.
✗ Branch 2 not taken.
 5936 cur_el.sons_.resize(4);
674 5936 number_of_elements_ += 3;
675 5936 max_level_ =
676
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5936 times.
 5936 max_level_ < cur_el.level_ + 1 ? cur_el.level_ + 1 : max_level_;
677 // auto it = leafs_.find(compute_global_id(cur_el));
678 // if (it != leafs_.end()) leafs_.erase(it);
679
2/2
✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
 29680 for (auto i = 0; i < 4; ++i) {
680 // add linked list structure to the panels
681 //////////////////////////////////////////////////////////////////////////
682
2/2
✓ Branch 0 taken 5936 times.
✓ Branch 1 taken 17808 times.
 23744 if (i == 0) {
683 5936 cur_el.sons_[i].prev_ = cur_el.prev_;
684
2/2
✓ Branch 0 taken 5695 times.
✓ Branch 1 taken 241 times.
 5936 if (cur_el.prev_ != nullptr)
685 5695 (cur_el.prev_)->next_ = std::addressof(cur_el.sons_[i]);
686 5936 cur_el.prev_ = nullptr;
687
2/2
✓ Branch 1 taken 241 times.
✓ Branch 2 taken 5695 times.
 5936 if (std::addressof(cur_el) == pfirst_)
688 241 pfirst_ = std::addressof(cur_el.sons_[i]);
689
690 } else {
691 17808 cur_el.sons_[i].prev_ = std::addressof(cur_el.sons_[i - 1]);
692 }
693
2/2
✓ Branch 0 taken 5936 times.
✓ Branch 1 taken 17808 times.
 23744 if (i == 3) {
694 5936 cur_el.sons_[i].next_ = cur_el.next_;
695
2/2
✓ Branch 0 taken 5695 times.
✓ Branch 1 taken 241 times.
 5936 if (cur_el.next_ != nullptr)
696 5695 (cur_el.next_)->prev_ = std::addressof(cur_el.sons_[i]);
697 5936 cur_el.next_ = nullptr;
698
2/2
✓ Branch 1 taken 241 times.
✓ Branch 2 taken 5695 times.
 5936 if (std::addressof(cur_el) == plast_)
699 241 plast_ = std::addressof(cur_el.sons_[i]);
700 } else {
701 17808 cur_el.sons_[i].next_ = std::addressof(cur_el.sons_[i + 1]);
702 }
703 //////////////////////////////////////////////////////////////////////////
704 23744 cur_el.sons_[i].patch_ = cur_el.patch_;
705 23744 cur_el.sons_[i].level_ = cur_el.level_ + 1;
706 23744 cur_el.sons_[i].id_ = 4 * cur_el.id_ + i;
707
1/2
✓ Branch 2 taken 23744 times.
✗ Branch 3 not taken.
 23744 cur_el.sons_[i].adjcents_ = std::vector<ElementTreeNode *>(4, nullptr);
708
1/2
✓ Branch 1 taken 23744 times.
✗ Branch 2 not taken.
 23744 cur_el.sons_[i].llc_(0) =
709
1/2
✓ Branch 1 taken 23744 times.
✗ Branch 2 not taken.
 23744 cur_el.llc_(0) + Constants::llcs[0][i] / double(1 << cur_el.level_);
710
1/2
✓ Branch 1 taken 23744 times.
✗ Branch 2 not taken.
 23744 cur_el.sons_[i].llc_(1) =
711
1/2
✓ Branch 1 taken 23744 times.
✗ Branch 2 not taken.
 23744 cur_el.llc_(1) + Constants::llcs[1][i] / double(1 << cur_el.level_);
712
1/2
✓ Branch 3 taken 23744 times.
✗ Branch 4 not taken.
 23744 elements.push_back(std::addressof(cur_el.sons_[i]));
713 }
714 // set vertices
715 // first element
716
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[0].vertices_.push_back(cur_el.vertices_[0]);
717
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[0].vertices_.push_back(ptIds[0]);
718
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[0].vertices_.push_back(ptIds[4]);
719
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[0].vertices_.push_back(ptIds[3]);
720 // second element
721
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[1].vertices_.push_back(ptIds[0]);
722
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[1].vertices_.push_back(cur_el.vertices_[1]);
723
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[1].vertices_.push_back(ptIds[1]);
724
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[1].vertices_.push_back(ptIds[4]);
725 // third element
726
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[2].vertices_.push_back(ptIds[4]);
727
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[2].vertices_.push_back(ptIds[1]);
728
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[2].vertices_.push_back(cur_el.vertices_[2]);
729
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[2].vertices_.push_back(ptIds[2]);
730 // fourth element
731
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[3].vertices_.push_back(ptIds[3]);
732
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[3].vertices_.push_back(ptIds[4]);
733
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[3].vertices_.push_back(ptIds[2]);
734
1/2
✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
 5936 cur_el.sons_[3].vertices_.push_back(cur_el.vertices_[3]);
735 // fix adjecency relations
736
1/2
✓ Branch 1 taken 5936 times.
✗ Branch 2 not taken.
 5936 updateTopology(elements);
737
738 5936 return;
739 5936 }
740 //////////////////////////////////////////////////////////////////////////////
741 /// member variables
742 //////////////////////////////////////////////////////////////////////////////
743 std::shared_ptr<PatchVector> geometry_;
744 ElementTreeNode root_;
745 ElementTreeNode *pfirst_;
746 ElementTreeNode *plast_;
747 int number_of_patches_;
748 int max_level_;
749 int number_of_points_;
750 int number_of_elements_;
751 //////////////////////////////////////////////////////////////////////////////
752 };
753 } // namespace Bembel
754 #endif // BEMBEL_SRC_CLUSTERTREE_ELEMENTTREE_HPP_
755