| 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_H2MATRIX_BLOCKCLUSTERTREE_HPP_ | ||
| 12 | #define BEMBEL_SRC_H2MATRIX_BLOCKCLUSTERTREE_HPP_ | ||
| 13 | |||
| 14 | namespace Bembel { | ||
| 15 | |||
| 16 | template <typename Scalar> | ||
| 17 | class BlockClusterTree; | ||
| 18 | |||
| 19 | struct BlockClusterAdmissibility { | ||
| 20 | enum { Refine = 0, LowRank = 1, Dense = 2 }; | ||
| 21 | }; | ||
| 22 | |||
| 23 | template <typename Scalar> | ||
| 24 | struct BlockClusterTreeParameters { | ||
| 25 | 5 | BlockClusterTreeParameters() | |
| 26 |
1/2✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
|
5 | : eta_(-1), min_cluster_level_(-1), max_level_(-1) {} |
| 27 | GaussSquare<Constants::maximum_quadrature_degree> GS_; | ||
| 28 | Eigen::Matrix<double, 12, Eigen::Dynamic> ffield_qnodes_; | ||
| 29 | double eta_; // eta from admissibility condition | ||
| 30 | int ffield_deg_; // todo @Michael comment this | ||
| 31 | int min_cluster_level_; // a tree leaf has 4^min_cluster_level_ elements | ||
| 32 | int max_level_; // depth of block cluster tree | ||
| 33 | int polynomial_degree_; // todo @Michael comment this | ||
| 34 | int polynomial_degree_plus_one_squared_; // todo @Michael comment this | ||
| 35 | const ElementTreeNode *et_root_; | ||
| 36 | }; | ||
| 37 | |||
| 38 | template <typename Scalar> | ||
| 39 | class BlockClusterTree { | ||
| 40 | public: | ||
| 41 | ////////////////////////////////////////////////////////////////////////////// | ||
| 42 | /// constructors | ||
| 43 | ////////////////////////////////////////////////////////////////////////////// | ||
| 44 | 3068 | BlockClusterTree() | |
| 45 | 3068 | : parameters_(nullptr), | |
| 46 | 3068 | sons_(GenericMatrix<BlockClusterTree>()), | |
| 47 | 3068 | leaf_(nullptr), | |
| 48 | 3068 | leaf_pointers_(nullptr), | |
| 49 | 3068 | cluster1_(nullptr), | |
| 50 | 3068 | cluster2_(nullptr), | |
| 51 | 3068 | rows_(0), | |
| 52 | 3068 | cols_(0), | |
| 53 | 3068 | cc_(-1) {} | |
| 54 | |||
| 55 | ✗ | BlockClusterTree(BlockClusterTree &&other) noexcept { | |
| 56 | ✗ | parameters_ = other.parameters_; | |
| 57 | ✗ | sons_ = std::move(other.sons_); | |
| 58 | ✗ | leaf_ = other.leaf_; | |
| 59 | ✗ | leaf_pointers_ = other.leaf_pointers_; | |
| 60 | ✗ | cluster1_ = other.cluster1_; | |
| 61 | ✗ | cluster2_ = other.cluster2_; | |
| 62 | ✗ | rows_ = other.rows_; | |
| 63 | ✗ | cols_ = other.cols_; | |
| 64 | ✗ | cc_ = other.cc_; | |
| 65 | ✗ | updateLeafPointers(); | |
| 66 | ✗ | } | |
| 67 | |||
| 68 | 4904 | BlockClusterTree(const BlockClusterTree &other) noexcept { | |
| 69 | 4904 | parameters_ = other.parameters_; | |
| 70 | 4904 | sons_ = other.sons_; | |
| 71 | // tree leaf is deep copied for const copy constructor! This has to be | ||
| 72 | // taken care of in some sense when going for Hmatrix arithmetics... | ||
| 73 |
2/2✓ Branch 1 taken 4608 times.
✓ Branch 2 taken 296 times.
|
4904 | if (other.leaf_) { |
| 74 | 4608 | leaf_ = std::make_shared< | |
| 75 | TreeLeaf<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>>(); | ||
| 76 | 4608 | *leaf_ = static_cast<const TreeLeaf< | |
| 77 | Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>> &>( | ||
| 78 | 4608 | *(other.leaf_)); | |
| 79 | } | ||
| 80 | 4904 | leaf_pointers_ = std::make_shared<std::vector<BlockClusterTree *>>(); | |
| 81 | 4904 | cluster1_ = other.cluster1_; | |
| 82 | 4904 | cluster2_ = other.cluster2_; | |
| 83 | 4904 | rows_ = other.rows_; | |
| 84 | 4904 | cols_ = other.cols_; | |
| 85 | 4904 | cc_ = other.cc_; | |
| 86 | 4904 | updateLeafPointers(); | |
| 87 | 4904 | } | |
| 88 | |||
| 89 | template <typename LinOp, typename AnsatzSpace> | ||
| 90 | 5 | BlockClusterTree(const LinOp &linear_operator, | |
| 91 | 5 | const AnsatzSpace &ansatz_space) { | |
| 92 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | init_BlockClusterTree(linear_operator, ansatz_space); |
| 93 | 5 | } | |
| 94 | ////////////////////////////////////////////////////////////////////////////// | ||
| 95 | /// assignment operator | ||
| 96 | ////////////////////////////////////////////////////////////////////////////// | ||
| 97 | 8 | BlockClusterTree &operator=(BlockClusterTree other) noexcept { | |
| 98 | 8 | std::swap(parameters_, other.parameters_); | |
| 99 | 8 | std::swap(sons_, other.sons_); | |
| 100 | 8 | std::swap(leaf_, other.leaf_); | |
| 101 | 8 | std::swap(leaf_pointers_, other.leaf_pointers_); | |
| 102 | 8 | std::swap(cluster1_, other.cluster1_); | |
| 103 | 8 | std::swap(cluster2_, other.cluster2_); | |
| 104 | 8 | std::swap(rows_, other.rows_); | |
| 105 | 8 | std::swap(cols_, other.cols_); | |
| 106 | 8 | std::swap(cc_, other.cc_); | |
| 107 | 8 | updateLeafPointers(); | |
| 108 | 8 | return *this; | |
| 109 | } | ||
| 110 | ////////////////////////////////////////////////////////////////////////////// | ||
| 111 | /// init | ||
| 112 | ////////////////////////////////////////////////////////////////////////////// | ||
| 113 | 5 | void set_parameters(double eta = 1.6, int min_cluster_level = 1) { | |
| 114 | 5 | parameters_ = std::make_shared<BlockClusterTreeParameters<Scalar>>(); | |
| 115 | 5 | parameters_->eta_ = eta; | |
| 116 | 5 | parameters_->min_cluster_level_ = min_cluster_level; | |
| 117 | 5 | return; | |
| 118 | } | ||
| 119 | |||
| 120 | template <typename LinOp, typename AnsatzSpace> | ||
| 121 | 5 | void init_BlockClusterTree(const LinOp &linOp, | |
| 122 | const AnsatzSpace &ansatz_space) { | ||
| 123 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | if (parameters_ == nullptr) set_parameters(); |
| 124 | // get element tree from ansatz space | ||
| 125 | const ElementTree &element_tree = | ||
| 126 | 5 | ansatz_space.get_superspace().get_mesh().get_element_tree(); | |
| 127 | 5 | parameters_->et_root_ = std::addressof(element_tree.root()); | |
| 128 | 5 | cluster1_ = std::addressof(element_tree.root()); | |
| 129 | 5 | cluster2_ = std::addressof(element_tree.root()); | |
| 130 | // set parameters for matrix assembly | ||
| 131 | 5 | parameters_->max_level_ = element_tree.get_max_level(); | |
| 132 | 5 | parameters_->polynomial_degree_ = | |
| 133 | 10 | ansatz_space.get_superspace().get_polynomial_degree(); | |
| 134 | 10 | parameters_->polynomial_degree_plus_one_squared_ = | |
| 135 | 5 | (parameters_->polynomial_degree_ + 1) * | |
| 136 | 5 | (parameters_->polynomial_degree_ + 1); | |
| 137 | 10 | parameters_->ffield_deg_ = | |
| 138 | 5 | linOp.get_FarfieldQuadratureDegree(parameters_->polynomial_degree_); | |
| 139 | // set up leaf_pointers | ||
| 140 | 5 | leaf_pointers_ = std::make_shared<std::vector<BlockClusterTree *>>(); | |
| 141 | 5 | leaf_pointers_->clear(); | |
| 142 | // block cluster tree assembly | ||
| 143 | 5 | rows_ = element_tree.get_number_of_elements(); | |
| 144 | 5 | cols_ = element_tree.get_number_of_elements(); | |
| 145 | // we let appendSubtree handle everything, since the root always | ||
| 146 | // returns 0 in compare cluster | ||
| 147 | 5 | appendSubtree(linOp, ansatz_space, cluster1_, cluster2_); | |
| 148 | 5 | updateLeafPointers(); | |
| 149 | 5 | return; | |
| 150 | } | ||
| 151 | ////////////////////////////////////////////////////////////////////////////// | ||
| 152 | /// methods | ||
| 153 | ////////////////////////////////////////////////////////////////////////////// | ||
| 154 | template <typename LinOp, typename AnsatzSpace> | ||
| 155 | 3065 | void appendSubtree(const LinOp &linear_operator, | |
| 156 | const AnsatzSpace &ansatz_space, | ||
| 157 | const ElementTreeNode *cluster1, | ||
| 158 | const ElementTreeNode *cluster2) { | ||
| 159 | 3065 | cc_ = compareCluster(*cluster1, *cluster2); | |
| 160 | // there are children to handle | ||
| 161 |
2/2✓ Branch 0 taken 185 times.
✓ Branch 1 taken 2880 times.
|
3065 | if (cc_ == BlockClusterAdmissibility::Refine) { |
| 162 | // reserve memory for sons_ | ||
| 163 | 185 | sons_.resize(cluster1->sons_.size(), cluster2->sons_.size()); | |
| 164 |
2/2✓ Branch 1 taken 750 times.
✓ Branch 2 taken 185 times.
|
935 | for (auto j = 0; j < sons_.cols(); ++j) |
| 165 |
2/2✓ Branch 1 taken 3060 times.
✓ Branch 2 taken 750 times.
|
3810 | for (auto i = 0; i < sons_.rows(); ++i) { |
| 166 | 3060 | const ElementTreeNode &son1 = cluster1->sons_[i]; | |
| 167 | 3060 | const ElementTreeNode &son2 = cluster2->sons_[j]; | |
| 168 | 3060 | sons_(i, j).parameters_ = parameters_; | |
| 169 | 3060 | sons_(i, j).cluster1_ = std::addressof(son1); | |
| 170 | 3060 | sons_(i, j).cluster2_ = std::addressof(son2); | |
| 171 | 3060 | sons_(i, j).rows_ = std::distance(son1.begin(), son1.end()); | |
| 172 | 3060 | sons_(i, j).cols_ = std::distance(son2.begin(), son2.end()); | |
| 173 | // let recursion handle the rest | ||
| 174 | 3060 | sons_(i, j).appendSubtree(linear_operator, ansatz_space, | |
| 175 | std::addressof(son1), std::addressof(son2)); | ||
| 176 | } | ||
| 177 | } else { | ||
| 178 | 2880 | leaf_ = std::make_shared< | |
| 179 | TreeLeaf<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>>(); | ||
| 180 | } | ||
| 181 | 3065 | return; | |
| 182 | } | ||
| 183 | |||
| 184 | void printNode() const { | ||
| 185 | std::cout << "{" << std::endl; | ||
| 186 | std::cout << "cluster: " << cluster1_ << " & " << cluster2_ << std::endl; | ||
| 187 | std::cout << "rows: " << rows_ << " cols: " << cols_ << std::endl; | ||
| 188 | std::cout << "format: " << cc_ << std::endl; | ||
| 189 | std::cout << "children: " << sons_.rows() << " x " << sons_.cols() | ||
| 190 | << std::endl; | ||
| 191 | std::cout << "leaf: " << leaf_ << std::endl; | ||
| 192 | std::cout << "parameters: " << parameters_ << std::endl; | ||
| 193 | |||
| 194 | std::cout << "}" << std::endl; | ||
| 195 | return; | ||
| 196 | } | ||
| 197 | ////////////////////////////////////////////////////////////////////////////// | ||
| 198 | /// getter | ||
| 199 | ////////////////////////////////////////////////////////////////////////////// | ||
| 200 | 448704 | int get_cc() const { return cc_; } | |
| 201 | int rows() const { return rows_; } | ||
| 202 | int cols() const { return cols_; } | ||
| 203 | 302976 | const ElementTreeNode *get_cluster1() const { return cluster1_; } | |
| 204 | 152928 | const ElementTreeNode *get_cluster2() const { return cluster2_; } | |
| 205 | int get_level() const { return get_cluster1()->get_level(); } | ||
| 206 | 8 | typename std::vector<BlockClusterTree *>::iterator lbegin() { | |
| 207 | 8 | return leaf_pointers_->begin(); | |
| 208 | } | ||
| 209 | 2885 | typename std::vector<BlockClusterTree *>::iterator lend() { | |
| 210 | 2885 | return leaf_pointers_->end(); | |
| 211 | } | ||
| 212 | 774 | typename std::vector<BlockClusterTree *>::const_iterator clbegin() const { | |
| 213 | 774 | return leaf_pointers_->begin(); | |
| 214 | } | ||
| 215 | 446598 | typename std::vector<BlockClusterTree *>::const_iterator clend() const { | |
| 216 | 446598 | return leaf_pointers_->end(); | |
| 217 | } | ||
| 218 | 591552 | int get_row_start_index() { | |
| 219 | 591552 | return get_parameters().polynomial_degree_plus_one_squared_ * | |
| 220 | 591552 | std::distance(get_parameters().et_root_->begin(), | |
| 221 | 1183104 | cluster1_->begin()); | |
| 222 | } | ||
| 223 | 295776 | int get_row_end_index() { | |
| 224 | 295776 | return get_parameters().polynomial_degree_plus_one_squared_ * | |
| 225 | 295776 | std::distance(get_parameters().et_root_->begin(), cluster1_->end()); | |
| 226 | } | ||
| 227 | 591552 | int get_col_start_index() { | |
| 228 | 591552 | return get_parameters().polynomial_degree_plus_one_squared_ * | |
| 229 | 591552 | std::distance(get_parameters().et_root_->begin(), | |
| 230 | 1183104 | cluster2_->begin()); | |
| 231 | } | ||
| 232 | 295776 | int get_col_end_index() { | |
| 233 | 295776 | return get_parameters().polynomial_degree_plus_one_squared_ * | |
| 234 | 295776 | std::distance(get_parameters().et_root_->begin(), cluster2_->end()); | |
| 235 | } | ||
| 236 | 300702 | const BlockClusterTreeParameters<Scalar> &get_parameters() const { | |
| 237 | 300702 | return *parameters_; | |
| 238 | } | ||
| 239 | 3549317 | BlockClusterTreeParameters<Scalar> &get_parameters() { return *parameters_; } | |
| 240 | |||
| 241 | const TreeLeaf<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>> | ||
| 242 | &get_leaf() const { | ||
| 243 | return *leaf_; | ||
| 244 | } | ||
| 245 | 452064 | TreeLeaf<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>> &get_leaf() { | |
| 246 | 452064 | return *leaf_; | |
| 247 | } | ||
| 248 | ////////////////////////////////////////////////////////////////////////////// | ||
| 249 | /// private members | ||
| 250 | ////////////////////////////////////////////////////////////////////////////// | ||
| 251 | private: | ||
| 252 | /** | ||
| 253 | * \brief determines admissible clusters | ||
| 254 | **/ | ||
| 255 | 3065 | int compareCluster(const ElementTreeNode &cluster1, | |
| 256 | const ElementTreeNode &cluster2) { | ||
| 257 | int r; | ||
| 258 |
1/2✓ Branch 2 taken 3065 times.
✗ Branch 3 not taken.
|
3065 | double dist = (cluster1.midpoint_ - cluster2.midpoint_).norm() - |
| 259 | 3065 | cluster1.radius_ - cluster2.radius_; | |
| 260 |
2/2✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 1145 times.
|
3065 | dist = dist > 0 ? dist : 0; |
| 261 |
2/2✓ Branch 0 taken 176 times.
✓ Branch 1 taken 2889 times.
|
3065 | double max_radius = cluster1.radius_ > cluster2.radius_ ? cluster1.radius_ |
| 262 | : cluster2.radius_; | ||
| 263 | |||
| 264 | /** | ||
| 265 | * \todo replace by 2*max_radius | ||
| 266 | */ | ||
| 267 |
5/6✓ Branch 0 taken 3065 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 2009 times.
✓ Branch 4 taken 1056 times.
✓ Branch 5 taken 2009 times.
✓ Branch 6 taken 1056 times.
|
3065 | if (dist < 0 || max_radius >= parameters_->eta_ * dist) |
| 268 | 2009 | r = BlockClusterAdmissibility::Refine; | |
| 269 | else | ||
| 270 | 1056 | r = BlockClusterAdmissibility::LowRank; | |
| 271 | |||
| 272 |
4/4✓ Branch 0 taken 2009 times.
✓ Branch 1 taken 1056 times.
✓ Branch 2 taken 1824 times.
✓ Branch 3 taken 1241 times.
|
5074 | if (r == BlockClusterAdmissibility::Refine && |
| 273 | 2009 | parameters_->max_level_ - cluster1.level_ <= | |
| 274 |
2/2✓ Branch 1 taken 1824 times.
✓ Branch 2 taken 185 times.
|
2009 | parameters_->min_cluster_level_) |
| 275 | 1824 | r = BlockClusterAdmissibility::Dense; | |
| 276 | |||
| 277 | 3065 | return r; | |
| 278 | } | ||
| 279 | 4917 | void updateLeafPointers() { | |
| 280 | // make sure we have the root of the block cluster tree calling this method | ||
| 281 |
3/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 4896 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 21 times.
|
4917 | if (cluster1_->id_ != -1 || cluster1_ != cluster2_) { |
| 282 | 4896 | return; | |
| 283 | } else { | ||
| 284 | 21 | leaf_pointers_->clear(); | |
| 285 | 21 | updateLeafPointers_recursion(*this); | |
| 286 | } | ||
| 287 | } | ||
| 288 | 12873 | void updateLeafPointers_recursion(BlockClusterTree &child) { | |
| 289 |
2/2✓ Branch 0 taken 777 times.
✓ Branch 1 taken 12096 times.
|
12873 | if (child.cc_ == BlockClusterAdmissibility::Refine) { |
| 290 |
2/2✓ Branch 1 taken 3150 times.
✓ Branch 2 taken 777 times.
|
3927 | for (auto j = 0; j < child.sons_.cols(); ++j) |
| 291 |
2/2✓ Branch 1 taken 12852 times.
✓ Branch 2 taken 3150 times.
|
16002 | for (auto i = 0; i < child.sons_.rows(); ++i) |
| 292 | 12852 | updateLeafPointers_recursion(child.sons_(i, j)); | |
| 293 | } else { | ||
| 294 |
1/2✓ Branch 3 taken 12096 times.
✗ Branch 4 not taken.
|
12096 | leaf_pointers_->push_back(std::addressof(child)); |
| 295 | } | ||
| 296 | 12873 | return; | |
| 297 | } | ||
| 298 | ////////////////////////////////////////////////////////////////////////////// | ||
| 299 | /// member variables | ||
| 300 | ////////////////////////////////////////////////////////////////////////////// | ||
| 301 | std::shared_ptr<BlockClusterTreeParameters<Scalar>> parameters_; | ||
| 302 | GenericMatrix<BlockClusterTree> sons_; | ||
| 303 | std::shared_ptr< | ||
| 304 | TreeLeaf<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>> | ||
| 305 | leaf_; | ||
| 306 | std::shared_ptr<std::vector<BlockClusterTree *>> leaf_pointers_; | ||
| 307 | const ElementTreeNode *cluster1_; | ||
| 308 | const ElementTreeNode *cluster2_; | ||
| 309 | int rows_; | ||
| 310 | int cols_; | ||
| 311 | int cc_; | ||
| 312 | }; | ||
| 313 | |||
| 314 | } // namespace Bembel | ||
| 315 | #endif // BEMBEL_SRC_H2MATRIX_BLOCKCLUSTERTREE_HPP_ | ||
| 316 |