| 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 725 times.
✓ Branch 1 taken 168 times.
|
893 | for (auto i = 0; i < number_of_patches_; ++i) { |
| 120 |
1/2✓ Branch 2 taken 725 times.
✗ Branch 3 not taken.
|
725 | root.sons_[i].adjcents_ = std::vector<ElementTreeNode *>(4, nullptr); |
| 121 |
1/2✓ Branch 2 taken 725 times.
✗ Branch 3 not taken.
|
725 | root.sons_[i].vertices_.resize(4); |
| 122 | 725 | root.sons_[i].id_ = i; | |
| 123 | 725 | root.sons_[i].level_ = 0; | |
| 124 | 725 | root.sons_[i].patch_ = i; | |
| 125 | // add linked list structure to the panels | ||
| 126 |
2/2✓ Branch 0 taken 168 times.
✓ Branch 1 taken 557 times.
|
725 | if (i == 0) { |
| 127 | 168 | root.sons_[0].prev_ = nullptr; | |
| 128 | 168 | pfirst_ = std::addressof(root.sons_[0]); | |
| 129 | } else { | ||
| 130 | 557 | root.sons_[i].prev_ = std::addressof(root.sons_[i - 1]); | |
| 131 | } | ||
| 132 |
2/2✓ Branch 0 taken 168 times.
✓ Branch 1 taken 557 times.
|
725 | if (i == number_of_patches_ - 1) { |
| 133 | 168 | root.sons_[i].next_ = nullptr; | |
| 134 | 168 | plast_ = std::addressof(root.sons_[i]); | |
| 135 | } else { | ||
| 136 | 557 | 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 2900 times.
✓ Branch 1 taken 725 times.
|
3625 | for (auto j = 0; j < 4; ++j) { |
| 141 | 2900 | v = (*geometry_)[i].eval(Constants::corners[0][j], | |
| 142 |
1/2✓ Branch 1 taken 2900 times.
✗ Branch 2 not taken.
|
2900 | Constants::corners[1][j]); |
| 143 | 2900 | unsigned int index = 0; | |
| 144 |
2/2✓ Branch 1 taken 12094 times.
✓ Branch 2 taken 1732 times.
|
13826 | for (; index < uniquePts.size(); ++index) |
| 145 |
4/6✓ Branch 2 taken 12094 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12094 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1168 times.
✓ Branch 8 taken 10926 times.
|
12094 | if ((uniquePts[index] - v).norm() < Constants::pt_comp_tolerance) |
| 146 | 1168 | break; | |
| 147 |
2/2✓ Branch 1 taken 1168 times.
✓ Branch 2 taken 1732 times.
|
2900 | if (index != uniquePts.size()) { |
| 148 | 1168 | root.sons_[i].vertices_[j] = index; | |
| 149 | } else { | ||
| 150 |
1/2✓ Branch 1 taken 1732 times.
✗ Branch 2 not taken.
|
1732 | uniquePts.push_back(v); |
| 151 | 1732 | root.sons_[i].vertices_[j] = number_of_points_; | |
| 152 | 1732 | ++number_of_points_; | |
| 153 | } | ||
| 154 | } | ||
| 155 |
1/2✓ Branch 3 taken 725 times.
✗ Branch 4 not taken.
|
725 | 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 18533 times.
✓ Branch 5 taken 168 times.
|
18701 | for (auto it = pbegin(); it != pend(); ++it) { |
| 213 | 18533 | double h = double(1) / double(1 << it->level_); | |
| 214 |
1/2✓ Branch 3 taken 18533 times.
✗ Branch 4 not taken.
|
18533 | pts.col(it->vertices_[0]) = |
| 215 |
4/8✓ Branch 5 taken 18533 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18533 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18533 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18533 times.
✗ Branch 16 not taken.
|
37066 | (*geometry_)[it->patch_].eval(it->llc_(0), it->llc_(1)); |
| 216 |
1/2✓ Branch 3 taken 18533 times.
✗ Branch 4 not taken.
|
18533 | pts.col(it->vertices_[1]) = |
| 217 |
4/8✓ Branch 5 taken 18533 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18533 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18533 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18533 times.
✗ Branch 16 not taken.
|
37066 | (*geometry_)[it->patch_].eval(it->llc_(0) + h, it->llc_(1)); |
| 218 |
1/2✓ Branch 3 taken 18533 times.
✗ Branch 4 not taken.
|
18533 | pts.col(it->vertices_[2]) = |
| 219 |
4/8✓ Branch 5 taken 18533 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18533 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18533 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18533 times.
✗ Branch 16 not taken.
|
37066 | (*geometry_)[it->patch_].eval(it->llc_(0) + h, it->llc_(1) + h); |
| 220 |
1/2✓ Branch 3 taken 18533 times.
✗ Branch 4 not taken.
|
18533 | pts.col(it->vertices_[3]) = |
| 221 |
4/8✓ Branch 5 taken 18533 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18533 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 18533 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 18533 times.
✗ Branch 16 not taken.
|
37066 | (*geometry_)[it->patch_].eval(it->llc_(0), it->llc_(1) + h); |
| 222 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18533 times.
|
18533 | 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 | 99825469 | 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 | 69718 | ElementTreeNode::const_iterator pend() const { | |
| 308 | 69718 | 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 | 36851 | 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 725 times.
✓ Branch 2 taken 168 times.
|
893 | for (auto i = 0; i < root_.sons_.size(); ++i) |
| 358 |
1/2✓ Branch 2 taken 725 times.
✗ Branch 3 not taken.
|
725 | 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 | 24469 | void computeElementEnclosings_recursion(ElementTreeNode &el, | |
| 369 | const Eigen::MatrixXd &P) { | ||
| 370 |
2/4✓ Branch 1 taken 24469 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24469 times.
✗ Branch 5 not taken.
|
24469 | Eigen::Vector3d mp1, mp2; |
| 371 | double r1, r2; | ||
| 372 | // compute point list | ||
| 373 |
2/2✓ Branch 1 taken 18533 times.
✓ Branch 2 taken 5936 times.
|
24469 | if (!el.sons_.size()) { |
| 374 | // assign enclosing balls to leafs | ||
| 375 |
4/8✓ Branch 1 taken 18533 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 18533 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 18533 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 18533 times.
✗ Branch 12 not taken.
|
18533 | computeEnclosingBall(&mp1, &r1, P.col(el.vertices_[0]), 0, |
| 376 |
1/2✓ Branch 2 taken 18533 times.
✗ Branch 3 not taken.
|
18533 | P.col(el.vertices_[2]), 0); |
| 377 |
4/8✓ Branch 1 taken 18533 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 18533 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 18533 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 18533 times.
✗ Branch 12 not taken.
|
18533 | computeEnclosingBall(&mp2, &r2, P.col(el.vertices_[1]), 0, |
| 378 |
1/2✓ Branch 2 taken 18533 times.
✗ Branch 3 not taken.
|
18533 | P.col(el.vertices_[3]), 0); |
| 379 |
1/2✓ Branch 1 taken 18533 times.
✗ Branch 2 not taken.
|
18533 | computeEnclosingBall(&(el.midpoint_), &(el.radius_), mp1, r1, mp2, r2); |
| 380 | } else { | ||
| 381 | // handle the four(!!!) children | ||
| 382 |
2/2✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
|
29680 | for (auto i = 0; i < 4; ++i) |
| 383 |
1/2✓ Branch 2 taken 23744 times.
✗ Branch 3 not taken.
|
23744 | computeElementEnclosings_recursion(el.sons_[i], P); |
| 384 | // assign enclosing balls to fathers bottom up | ||
| 385 |
1/2✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
|
5936 | computeEnclosingBall(&mp1, &r1, el.sons_[0].midpoint_, |
| 386 | 5936 | el.sons_[0].radius_, el.sons_[2].midpoint_, | |
| 387 | 5936 | el.sons_[2].radius_); | |
| 388 |
1/2✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
|
5936 | computeEnclosingBall(&mp2, &r2, el.sons_[1].midpoint_, |
| 389 | 5936 | el.sons_[1].radius_, el.sons_[3].midpoint_, | |
| 390 | 5936 | el.sons_[3].radius_); | |
| 391 |
1/2✓ Branch 1 taken 5936 times.
✗ Branch 2 not taken.
|
5936 | computeEnclosingBall(&(el.midpoint_), &(el.radius_), mp1, r1, mp2, r2); |
| 392 | } | ||
| 393 | 48938 | return; | |
| 394 | } | ||
| 395 | ////////////////////////////////////////////////////////////////////////////// | ||
| 396 | /** | ||
| 397 | * \brief Prints all Elements of the Tree. | ||
| 398 | */ | ||
| 399 | void printPanels() const { | ||
| 400 | auto i = 0; | ||
| 401 | for (auto it = pbegin(); it != pend(); ++it) { | ||
| 402 | std::cout << "element[" << i++ << "] = "; | ||
| 403 | it->print(); | ||
| 404 | } | ||
| 405 | return; | ||
| 406 | } | ||
| 407 | ////////////////////////////////////////////////////////////////////////////// | ||
| 408 | /// other Stuff | ||
| 409 | ////////////////////////////////////////////////////////////////////////////// | ||
| 410 | /** | ||
| 411 | * \brief Return a matrix with all midpoints of the elements. | ||
| 412 | * | ||
| 413 | * \return The midpoint list as 3xN matrix with N the number of points. | ||
| 414 | */ | ||
| 415 | Eigen::MatrixXd generateMidpointList() const { | ||
| 416 | Eigen::MatrixXd retval(3, number_of_elements_); | ||
| 417 | unsigned int i = 0; | ||
| 418 | for (auto it = pbegin(); it != pend(); ++it) { | ||
| 419 | retval.col(i) = it->midpoint_; | ||
| 420 | ++i; | ||
| 421 | } | ||
| 422 | return retval; | ||
| 423 | } | ||
| 424 | ////////////////////////////////////////////////////////////////////////////// | ||
| 425 | /** | ||
| 426 | * \brief Return a matrix with all radii of the element enclosing. | ||
| 427 | * | ||
| 428 | * \return A vector containing the radii of the element enclosing. | ||
| 429 | */ | ||
| 430 | Eigen::MatrixXd generateRadiusList() const { | ||
| 431 | Eigen::VectorXd retval(number_of_elements_); | ||
| 432 | unsigned int i = 0; | ||
| 433 | for (auto it = pbegin(); it != pend(); ++it) { | ||
| 434 | retval(i) = it->radius_; | ||
| 435 | ++i; | ||
| 436 | } | ||
| 437 | return retval; | ||
| 438 | } | ||
| 439 | ////////////////////////////////////////////////////////////////////////////// | ||
| 440 | /** | ||
| 441 | * \brief Return a vector with computed global indices of the elements. | ||
| 442 | * | ||
| 443 | * \return A vector with global indices of the elements. | ||
| 444 | */ | ||
| 445 | Eigen::VectorXi generateElementLabels() const { | ||
| 446 | Eigen::VectorXi retval(number_of_elements_); | ||
| 447 | unsigned int i = 0; | ||
| 448 | for (auto it = pbegin(); it != pend(); ++it) { | ||
| 449 | retval(i) = compute_global_id(*it); | ||
| 450 | ++i; | ||
| 451 | } | ||
| 452 | return retval; | ||
| 453 | } | ||
| 454 | ////////////////////////////////////////////////////////////////////////////// | ||
| 455 | /** | ||
| 456 | * \brief Generate list of labels if elements are on the patch boundary or at | ||
| 457 | * the boundary of the geometry. | ||
| 458 | * | ||
| 459 | * The value at the index of the element is 1 if the element is at the patch | ||
| 460 | * boundary and there is a neighbor patch. If there is no neighbor patch then | ||
| 461 | * the value is -1. | ||
| 462 | * | ||
| 463 | * \return A vector of integers of size N with N the number of elements. | ||
| 464 | */ | ||
| 465 | Eigen::VectorXi generatePatchBoundaryLabels() const { | ||
| 466 | Eigen::VectorXi retval(number_of_elements_); | ||
| 467 | retval.setZero(); | ||
| 468 | unsigned int i = 0; | ||
| 469 | for (auto it = pbegin(); it != pend(); ++it) { | ||
| 470 | for (auto j = 0; j < 4; ++j) | ||
| 471 | if (it->adjcents_[j] == nullptr) { | ||
| 472 | retval[i] = -1; | ||
| 473 | break; | ||
| 474 | } else if (it->adjcents_[j]->patch_ != it->patch_) { | ||
| 475 | ++(retval[i]); | ||
| 476 | } | ||
| 477 | ++i; | ||
| 478 | } | ||
| 479 | return retval; | ||
| 480 | } | ||
| 481 | ////////////////////////////////////////////////////////////////////////////// | ||
| 482 | /** | ||
| 483 | * \brief Generate list of labels if elements contained within a given patch. | ||
| 484 | * | ||
| 485 | * The value at the index of the element is 1 if the element is contained | ||
| 486 | * within the given patch. Otherwise its zero. | ||
| 487 | * | ||
| 488 | * \return A vector of integers of size N with N the number of elements. | ||
| 489 | */ | ||
| 490 | Eigen::VectorXi identifyPatch(unsigned int pn) const { | ||
| 491 | Eigen::VectorXi retval(number_of_elements_); | ||
| 492 | retval.setZero(); | ||
| 493 | assert(pn < number_of_patches_); | ||
| 494 | unsigned int i = 0; | ||
| 495 | for (auto it = pbegin(); it != pend(); ++it) { | ||
| 496 | retval[i] = it->patch_ == pn ? 1 : 0; | ||
| 497 | ++i; | ||
| 498 | } | ||
| 499 | return retval; | ||
| 500 | } | ||
| 501 | ////////////////////////////////////////////////////////////////////////////// | ||
| 502 | /// static members | ||
| 503 | ////////////////////////////////////////////////////////////////////////////// | ||
| 504 | /** | ||
| 505 | * \brief computes a ball enclosing the union of \f$B_r1(mp1)\f$ and \f$B_r2(mp2)\f$, | ||
| 506 | * i.e \f$B(mp,r)\supset B_r1(mp1) \cup B_r2(mp2)\f$. | ||
| 507 | */ | ||
| 508 | 73407 | static void computeEnclosingBall(Eigen::Vector3d *mp, double *r, | |
| 509 | const Eigen::Vector3d &mp1, double r1, | ||
| 510 | const Eigen::Vector3d &mp2, double r2) { | ||
| 511 | // compute distance vector of the two spheres | ||
| 512 |
1/2✓ Branch 1 taken 73407 times.
✗ Branch 2 not taken.
|
73407 | auto z = mp1 - mp2; |
| 513 |
2/4✓ Branch 1 taken 73407 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 73407 times.
✗ Branch 5 not taken.
|
73407 | auto norm = (mp1 - mp2).norm(); |
| 514 | // B(d2,r2) subset B(d1,r1) | ||
| 515 |
2/2✓ Branch 0 taken 24021 times.
✓ Branch 1 taken 49386 times.
|
73407 | if (norm + r2 <= r1) { |
| 516 |
1/2✓ Branch 1 taken 24021 times.
✗ Branch 2 not taken.
|
24021 | *mp = mp1; |
| 517 | 24021 | *r = r1; | |
| 518 | // B(d1,r1) subset B(d2,r2) | ||
| 519 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 49338 times.
|
49386 | } else if (norm + r1 <= r2) { |
| 520 |
1/2✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
|
48 | *mp = mp2; |
| 521 | 48 | *r = r2; | |
| 522 | // the union is not a ball | ||
| 523 | } else { | ||
| 524 |
5/10✓ Branch 1 taken 49338 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49338 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 49338 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 49338 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 49338 times.
✗ Branch 14 not taken.
|
49338 | *mp = 0.5 * (mp1 + mp2 + (r1 - r2) / norm * z); |
| 525 | 49338 | *r = 0.5 * (r1 + r2 + norm); | |
| 526 | 49338 | *r = 0.5 * (r1 + r2 + norm); | |
| 527 | } | ||
| 528 | 146814 | return; | |
| 529 | } | ||
| 530 | /** | ||
| 531 | * \brief Resolves neighborhood relations of the patches. | ||
| 532 | * | ||
| 533 | * \return A vector where each entry defines a patch interface or boundary. | ||
| 534 | * The entries correspond to [patchIndex1, edgeCase1, patchIndex2, edgeCase2]. | ||
| 535 | */ | ||
| 536 | ////////////////////////////////////////////////////////////////////////////// | ||
| 537 | 312 | std::vector<std::array<int, 4>> patchTopologyInfo() const { | |
| 538 | 312 | std::vector<std::array<int, 4>> retval; | |
| 539 |
2/2✓ Branch 4 taken 1426 times.
✓ Branch 5 taken 312 times.
|
1738 | for (auto it = root_.sons_.begin(); it != root_.sons_.end(); ++it) { |
| 540 |
2/2✓ Branch 0 taken 5704 times.
✓ Branch 1 taken 1426 times.
|
7130 | for (auto j = 0; j < 4; ++j) { |
| 541 | // do we have a neighbour? | ||
| 542 |
2/2✓ Branch 2 taken 2480 times.
✓ Branch 3 taken 3224 times.
|
5704 | if (it->adjcents_[j] != nullptr) { |
| 543 | 2480 | const ElementTreeNode &cur_neighbour = *(it->adjcents_[j]); | |
| 544 | // add the edge only if it->id_ < neighbour->id_ | ||
| 545 | // otherwise the edge has already been added | ||
| 546 |
2/2✓ Branch 1 taken 1240 times.
✓ Branch 2 taken 1240 times.
|
2480 | if (it->id_ < cur_neighbour.id_) { |
| 547 | 1240 | int k = 0; | |
| 548 |
1/2✓ Branch 0 taken 3086 times.
✗ Branch 1 not taken.
|
3086 | for (; k < 4; ++k) |
| 549 |
4/4✓ Branch 1 taken 1550 times.
✓ Branch 2 taken 1536 times.
✓ Branch 3 taken 1240 times.
✓ Branch 4 taken 1846 times.
|
4636 | if (cur_neighbour.adjcents_[k] != nullptr && |
| 550 |
2/2✓ Branch 2 taken 1240 times.
✓ Branch 3 taken 310 times.
|
1550 | cur_neighbour.adjcents_[k]->id_ == it->id_) |
| 551 | 1240 | break; | |
| 552 |
1/2✓ Branch 2 taken 1240 times.
✗ Branch 3 not taken.
|
1240 | retval.push_back({it->id_, cur_neighbour.id_, j, k}); |
| 553 | } | ||
| 554 | } else { | ||
| 555 |
1/2✓ Branch 2 taken 3224 times.
✗ Branch 3 not taken.
|
3224 | retval.push_back({it->id_, -1, j, -1}); |
| 556 | } | ||
| 557 | } | ||
| 558 | } | ||
| 559 | 312 | return retval; | |
| 560 | ✗ | } | |
| 561 | /** | ||
| 562 | * \brief The ordering of elements in the element tree does not correspond to | ||
| 563 | * the element order underlying the coefficient vector. This reordering can be | ||
| 564 | * computed for look ups by this function. | ||
| 565 | * | ||
| 566 | * Limitation to the uniform case! | ||
| 567 | * | ||
| 568 | * \return Vector with the tensor product index of the elements. | ||
| 569 | */ | ||
| 570 | ////////////////////////////////////////////////////////////////////////////// | ||
| 571 | 138 | std::vector<int> computeReorderingVector() const { | |
| 572 |
1/2✓ Branch 2 taken 138 times.
✗ Branch 3 not taken.
|
138 | std::vector<int> out(number_of_elements_); |
| 573 |
2/2✓ Branch 3 taken 14028 times.
✓ Branch 4 taken 138 times.
|
14166 | for (auto it = pbegin(); it != pend(); ++it) { |
| 574 | 14028 | const double h = it->get_h(); | |
| 575 | const Eigen::Vector2d ref_mid = | ||
| 576 |
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); |
| 577 |
1/2✓ Branch 1 taken 14028 times.
✗ Branch 2 not taken.
|
14028 | const int x_idx = std::floor(ref_mid(0) / h); |
| 578 |
1/2✓ Branch 1 taken 14028 times.
✗ Branch 2 not taken.
|
14028 | const int y_idx = std::floor(ref_mid(1) / h); |
| 579 | 14028 | const int num_in_one_dir = 1 << max_level_; | |
| 580 | 14028 | const int num_per_patch = num_in_one_dir * num_in_one_dir; | |
| 581 | const int tp_idx = | ||
| 582 | 14028 | it->patch_ * num_per_patch + y_idx * num_in_one_dir + x_idx; | |
| 583 | 14028 | out[tp_idx] = it->id_; | |
| 584 | } | ||
| 585 | 138 | return out; | |
| 586 | ✗ | } | |
| 587 | ////////////////////////////////////////////////////////////////////////////// | ||
| 588 | size_t compute_global_id(const ElementTreeNode &el) const { | ||
| 589 | return number_of_patches_ * | ||
| 590 | (((1 << el.level_) * (1 << el.level_) - 1) / 3) + | ||
| 591 | el.id_; | ||
| 592 | } | ||
| 593 | ////////////////////////////////////////////////////////////////////////////// | ||
| 594 | /// private members | ||
| 595 | ////////////////////////////////////////////////////////////////////////////// | ||
| 596 | private: | ||
| 597 | ////////////////////////////////////////////////////////////////////////////// | ||
| 598 | /// methods | ||
| 599 | ////////////////////////////////////////////////////////////////////////////// | ||
| 600 | template <typename T> | ||
| 601 | struct isEqual { | ||
| 602 | bool operator()(const T &v1, const T &v2) const { | ||
| 603 | return ((v1 - v2).norm() < Constants::pt_comp_tolerance); | ||
| 604 | } | ||
| 605 | }; | ||
| 606 | /** | ||
| 607 | * \brief function to set up the local topology, i.e. the adjacents | ||
| 608 | * of a refined element | ||
| 609 | */ | ||
| 610 | ////////////////////////////////////////////////////////////////////////////// | ||
| 611 | 6104 | void updateTopology(const std::vector<ElementTreeNode *> &elements) { | |
| 612 | 6104 | std::map<std::array<int, 2>, ElementTreeNode *> edges; | |
| 613 | std::array<int, 2> e1, e2; | ||
| 614 | // generate edge list for all elements in question | ||
| 615 |
2/2✓ Branch 1 taken 43541 times.
✓ Branch 2 taken 6104 times.
|
49645 | for (auto i = 0; i < elements.size(); ++i) |
| 616 |
2/2✓ Branch 0 taken 174164 times.
✓ Branch 1 taken 43541 times.
|
217705 | for (auto j = 0; j < 4; ++j) { |
| 617 | // compute a unique id for each edge | ||
| 618 | 174164 | const int v1 = elements[i]->vertices_[j]; | |
| 619 | 174164 | const int v2 = elements[i]->vertices_[(j + 1) % 4]; | |
| 620 |
4/4✓ Branch 0 taken 87241 times.
✓ Branch 1 taken 86923 times.
✓ Branch 2 taken 87241 times.
✓ Branch 3 taken 86923 times.
|
174164 | e1 = {v1 < v2 ? v1 : v2, v1 < v2 ? v2 : v1}; |
| 621 | // perferm a look up if the edge is already existing. | ||
| 622 |
1/2✓ Branch 1 taken 174164 times.
✗ Branch 2 not taken.
|
174164 | auto it = edges.find(e1); |
| 623 | // if so, identify the two neighbours | ||
| 624 |
2/2✓ Branch 2 taken 53060 times.
✓ Branch 3 taken 121104 times.
|
174164 | if (it != edges.end()) { |
| 625 | 53060 | elements[i]->adjcents_[j] = it->second; | |
| 626 | // now find the edge also in the patch that added it to the set | ||
| 627 |
1/2✓ Branch 0 taken 156501 times.
✗ Branch 1 not taken.
|
156501 | for (auto k = 0; k < 4; ++k) { |
| 628 | 156501 | const int v3 = it->second->vertices_[k]; | |
| 629 | 156501 | const int v4 = it->second->vertices_[(k + 1) % 4]; | |
| 630 |
4/4✓ Branch 0 taken 86670 times.
✓ Branch 1 taken 69831 times.
✓ Branch 2 taken 86670 times.
✓ Branch 3 taken 69831 times.
|
156501 | e2 = {v3 < v4 ? v3 : v4, v3 < v4 ? v4 : v3}; |
| 631 |
3/4✓ Branch 1 taken 156501 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 53060 times.
✓ Branch 4 taken 103441 times.
|
156501 | if (e1 == e2) { |
| 632 | 53060 | it->second->adjcents_[k] = elements[i]; | |
| 633 | 53060 | break; | |
| 634 | } | ||
| 635 | } | ||
| 636 | // otherwise add the edge to the list | ||
| 637 | } else { | ||
| 638 |
2/4✓ Branch 2 taken 121104 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 121104 times.
✗ Branch 6 not taken.
|
121104 | edges.insert(std::make_pair(e1, elements[i])); |
| 639 | } | ||
| 640 | } | ||
| 641 | 12208 | return; | |
| 642 | 6104 | } | |
| 643 | ////////////////////////////////////////////////////////////////////////////// | ||
| 644 | /** | ||
| 645 | * \brief This function refines the given ElementTreeNode. | ||
| 646 | * | ||
| 647 | * The function introduces 4 new elements and takes care of the newly | ||
| 648 | * introduces vertices. Furthermore, all the neighborhood relation of all | ||
| 649 | * surrounding elements are resolved. | ||
| 650 | */ | ||
| 651 | 5936 | void refineLeaf(ElementTreeNode &cur_el) { | |
| 652 | // check if we have actually a panel | ||
| 653 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5936 times.
|
5936 | if (cur_el.sons_.size()) return; |
| 654 |
1/2✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
|
5936 | std::vector<int> ptIds(5); |
| 655 |
1/2✓ Branch 2 taken 5936 times.
✗ Branch 3 not taken.
|
5936 | std::vector<int> refNeighbours(4); |
| 656 | 5936 | std::vector<ElementTreeNode *> elements; | |
| 657 | // determine position of p with respect to its neighbours | ||
| 658 |
2/2✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
|
29680 | for (auto i = 0; i < 4; ++i) { |
| 659 | // is there a neighbour? | ||
| 660 |
2/2✓ Branch 1 taken 19072 times.
✓ Branch 2 taken 4672 times.
|
23744 | if (cur_el.adjcents_[i] != nullptr) { |
| 661 |
1/2✓ Branch 0 taken 47680 times.
✗ Branch 1 not taken.
|
47680 | for (auto j = 0; j < 4; ++j) |
| 662 |
2/2✓ Branch 3 taken 19072 times.
✓ Branch 4 taken 28608 times.
|
47680 | if (cur_el.adjcents_[i]->adjcents_[j] == std::addressof(cur_el)) { |
| 663 | 19072 | refNeighbours[i] = j; | |
| 664 | 19072 | break; | |
| 665 | } | ||
| 666 | } else { | ||
| 667 | 4672 | refNeighbours[i] = -1; | |
| 668 | } | ||
| 669 | } | ||
| 670 | // determine new points | ||
| 671 |
2/2✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
|
29680 | for (auto i = 0; i < 4; ++i) { |
| 672 | // is there a neighbour? | ||
| 673 |
2/2✓ Branch 1 taken 19072 times.
✓ Branch 2 taken 4672 times.
|
23744 | if (cur_el.adjcents_[i] != nullptr) { |
| 674 | 19072 | ElementTreeNode &ref_cur_neighbour = *(cur_el.adjcents_[i]); | |
| 675 | // is the neighbour already refined? | ||
| 676 |
2/2✓ Branch 1 taken 9536 times.
✓ Branch 2 taken 9536 times.
|
19072 | if (ref_cur_neighbour.sons_.size()) { |
| 677 | // this is the midpoint of the shared edge | ||
| 678 | 19072 | ptIds[i] = ref_cur_neighbour.sons_[refNeighbours[i]] | |
| 679 | 19072 | .vertices_[(refNeighbours[i] + 1) % 4]; | |
| 680 | // these are the two elements adjacent to the edge | ||
| 681 | 9536 | elements.push_back( | |
| 682 |
1/2✓ Branch 4 taken 9536 times.
✗ Branch 5 not taken.
|
9536 | std::addressof(ref_cur_neighbour.sons_[refNeighbours[i]])); |
| 683 |
1/2✓ Branch 2 taken 9536 times.
✗ Branch 3 not taken.
|
9536 | elements.push_back(std::addressof( |
| 684 | 9536 | ref_cur_neighbour.sons_[(refNeighbours[i] + 1) % 4])); | |
| 685 | } else { | ||
| 686 | // otherwise add the point id to the tree | ||
| 687 | 9536 | ptIds[i] = number_of_points_; | |
| 688 | 9536 | ++number_of_points_; | |
| 689 | } | ||
| 690 | } else { | ||
| 691 | // otherwise add the point id to the tree | ||
| 692 | 4672 | ptIds[i] = number_of_points_; | |
| 693 | 4672 | ++number_of_points_; | |
| 694 | } | ||
| 695 | } | ||
| 696 | // add midpoint of the current element, which is always a new point | ||
| 697 | 5936 | ptIds[4] = number_of_points_; | |
| 698 | 5936 | ++number_of_points_; | |
| 699 | // set up new elements | ||
| 700 |
1/2✓ Branch 1 taken 5936 times.
✗ Branch 2 not taken.
|
5936 | cur_el.sons_.resize(4); |
| 701 | 5936 | number_of_elements_ += 3; | |
| 702 | 5936 | max_level_ = | |
| 703 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5936 times.
|
5936 | max_level_ < cur_el.level_ + 1 ? cur_el.level_ + 1 : max_level_; |
| 704 | // auto it = leafs_.find(compute_global_id(cur_el)); | ||
| 705 | // if (it != leafs_.end()) leafs_.erase(it); | ||
| 706 |
2/2✓ Branch 0 taken 23744 times.
✓ Branch 1 taken 5936 times.
|
29680 | for (auto i = 0; i < 4; ++i) { |
| 707 | // add linked list structure to the panels | ||
| 708 | ////////////////////////////////////////////////////////////////////////// | ||
| 709 |
2/2✓ Branch 0 taken 5936 times.
✓ Branch 1 taken 17808 times.
|
23744 | if (i == 0) { |
| 710 | 5936 | cur_el.sons_[i].prev_ = cur_el.prev_; | |
| 711 |
2/2✓ Branch 0 taken 5695 times.
✓ Branch 1 taken 241 times.
|
5936 | if (cur_el.prev_ != nullptr) |
| 712 | 5695 | (cur_el.prev_)->next_ = std::addressof(cur_el.sons_[i]); | |
| 713 | 5936 | cur_el.prev_ = nullptr; | |
| 714 |
2/2✓ Branch 1 taken 241 times.
✓ Branch 2 taken 5695 times.
|
5936 | if (std::addressof(cur_el) == pfirst_) |
| 715 | 241 | pfirst_ = std::addressof(cur_el.sons_[i]); | |
| 716 | |||
| 717 | } else { | ||
| 718 | 17808 | cur_el.sons_[i].prev_ = std::addressof(cur_el.sons_[i - 1]); | |
| 719 | } | ||
| 720 |
2/2✓ Branch 0 taken 5936 times.
✓ Branch 1 taken 17808 times.
|
23744 | if (i == 3) { |
| 721 | 5936 | cur_el.sons_[i].next_ = cur_el.next_; | |
| 722 |
2/2✓ Branch 0 taken 5695 times.
✓ Branch 1 taken 241 times.
|
5936 | if (cur_el.next_ != nullptr) |
| 723 | 5695 | (cur_el.next_)->prev_ = std::addressof(cur_el.sons_[i]); | |
| 724 | 5936 | cur_el.next_ = nullptr; | |
| 725 |
2/2✓ Branch 1 taken 241 times.
✓ Branch 2 taken 5695 times.
|
5936 | if (std::addressof(cur_el) == plast_) |
| 726 | 241 | plast_ = std::addressof(cur_el.sons_[i]); | |
| 727 | } else { | ||
| 728 | 17808 | cur_el.sons_[i].next_ = std::addressof(cur_el.sons_[i + 1]); | |
| 729 | } | ||
| 730 | ////////////////////////////////////////////////////////////////////////// | ||
| 731 | 23744 | cur_el.sons_[i].patch_ = cur_el.patch_; | |
| 732 | 23744 | cur_el.sons_[i].level_ = cur_el.level_ + 1; | |
| 733 | 23744 | cur_el.sons_[i].id_ = 4 * cur_el.id_ + i; | |
| 734 |
1/2✓ Branch 2 taken 23744 times.
✗ Branch 3 not taken.
|
23744 | cur_el.sons_[i].adjcents_ = std::vector<ElementTreeNode *>(4, nullptr); |
| 735 |
1/2✓ Branch 1 taken 23744 times.
✗ Branch 2 not taken.
|
23744 | cur_el.sons_[i].llc_(0) = |
| 736 |
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_); |
| 737 |
1/2✓ Branch 1 taken 23744 times.
✗ Branch 2 not taken.
|
23744 | cur_el.sons_[i].llc_(1) = |
| 738 |
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_); |
| 739 |
1/2✓ Branch 3 taken 23744 times.
✗ Branch 4 not taken.
|
23744 | elements.push_back(std::addressof(cur_el.sons_[i])); |
| 740 | } | ||
| 741 | // set vertices | ||
| 742 | // first element | ||
| 743 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[0].vertices_.push_back(cur_el.vertices_[0]); |
| 744 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[0].vertices_.push_back(ptIds[0]); |
| 745 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[0].vertices_.push_back(ptIds[4]); |
| 746 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[0].vertices_.push_back(ptIds[3]); |
| 747 | // second element | ||
| 748 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[1].vertices_.push_back(ptIds[0]); |
| 749 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[1].vertices_.push_back(cur_el.vertices_[1]); |
| 750 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[1].vertices_.push_back(ptIds[1]); |
| 751 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[1].vertices_.push_back(ptIds[4]); |
| 752 | // third element | ||
| 753 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[2].vertices_.push_back(ptIds[4]); |
| 754 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[2].vertices_.push_back(ptIds[1]); |
| 755 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[2].vertices_.push_back(cur_el.vertices_[2]); |
| 756 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[2].vertices_.push_back(ptIds[2]); |
| 757 | // fourth element | ||
| 758 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[3].vertices_.push_back(ptIds[3]); |
| 759 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[3].vertices_.push_back(ptIds[4]); |
| 760 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[3].vertices_.push_back(ptIds[2]); |
| 761 |
1/2✓ Branch 3 taken 5936 times.
✗ Branch 4 not taken.
|
5936 | cur_el.sons_[3].vertices_.push_back(cur_el.vertices_[3]); |
| 762 | // fix adjecency relations | ||
| 763 |
1/2✓ Branch 1 taken 5936 times.
✗ Branch 2 not taken.
|
5936 | updateTopology(elements); |
| 764 | |||
| 765 | 5936 | return; | |
| 766 | 5936 | } | |
| 767 | ////////////////////////////////////////////////////////////////////////////// | ||
| 768 | /// member variables | ||
| 769 | ////////////////////////////////////////////////////////////////////////////// | ||
| 770 | std::shared_ptr<PatchVector> geometry_; | ||
| 771 | ElementTreeNode root_; | ||
| 772 | ElementTreeNode *pfirst_; | ||
| 773 | ElementTreeNode *plast_; | ||
| 774 | int number_of_patches_; | ||
| 775 | int max_level_; | ||
| 776 | int number_of_points_; | ||
| 777 | int number_of_elements_; | ||
| 778 | ////////////////////////////////////////////////////////////////////////////// | ||
| 779 | }; | ||
| 780 | } // namespace Bembel | ||
| 781 | #endif // BEMBEL_SRC_CLUSTERTREE_ELEMENTTREE_HPP_ | ||
| 782 |