| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // This file is part of Bembel, the higher order C++ boundary element library. | ||
| 2 | // | ||
| 3 | // Copyright (C) 2022 see <http://www.bembel.eu> | ||
| 4 | // | ||
| 5 | // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, | ||
| 6 | // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, | ||
| 7 | // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This | ||
| 8 | // source code is subject to the GNU General Public License version 3 and | ||
| 9 | // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further | ||
| 10 | // information. | ||
| 11 | #ifndef BEMBEL_SRC_ANSATZSPACE_GLUE_HPP_ | ||
| 12 | #define BEMBEL_SRC_ANSATZSPACE_GLUE_HPP_ | ||
| 13 | |||
| 14 | namespace Bembel { | ||
| 15 | /* | ||
| 16 | |||
| 17 | Every edge has an index between 0 and 3, where | ||
| 18 | 0 = (0,0)->(1,0) | ||
| 19 | 1 = (1,0)->(1,1) | ||
| 20 | 2 = (1,1)->(0,1) | ||
| 21 | 3 = (0,1)->(0,0) | ||
| 22 | */ | ||
| 23 | |||
| 24 | /* | ||
| 25 | Example of Dof enumeration on one patch: | ||
| 26 | dimXdir = 4; | ||
| 27 | dimYdir = 3; | ||
| 28 | Edge 2 (1,1) | ||
| 29 | |||
| 30 | y |8 9 10 11| | ||
| 31 | A Edge 3 |4 5 6 7 | Edge 1 | ||
| 32 | | |0 1 2 3 | | ||
| 33 | |||
| 34 | Edge 0 | ||
| 35 | (0,0) -> x | ||
| 36 | */ | ||
| 37 | |||
| 38 | namespace GlueRoutines { | ||
| 39 | /** | ||
| 40 | * \brief Helper struct. | ||
| 41 | */ | ||
| 42 | struct dofIdentification { | ||
| 43 | std::vector<int> dofs; | ||
| 44 | int coef; | ||
| 45 | }; | ||
| 46 | /** | ||
| 47 | * \brief Helper struct for assembling the Glue for the different cases. | ||
| 48 | */ | ||
| 49 | template <typename Derived, unsigned int DF> | ||
| 50 | struct glue_identificationmaker_ { | ||
| 51 | static std::vector<dofIdentification> makeIdentification( | ||
| 52 | const std::vector<std::array<int, 4>> &edges_, | ||
| 53 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 54 | assert(false && "Needs to be specialized"); | ||
| 55 | return {}; | ||
| 56 | } | ||
| 57 | }; | ||
| 58 | |||
| 59 | } // namespace GlueRoutines | ||
| 60 | |||
| 61 | /** | ||
| 62 | * \ingroup AnsatzSpace | ||
| 63 | * \brief This class takes care of identifying DOFs on different edges, which | ||
| 64 | *must be identified with one another. | ||
| 65 | * | ||
| 66 | * The Glue class has routines to assemble continuous B-splines across patch | ||
| 67 | *boundaries depending on the template parameter | ||
| 68 | *LinearOperatorTraits<Derived>::Form. If this template parameter is equal to | ||
| 69 | *DifferentialForm::Continuous, the B-splines are globally glued continuously. | ||
| 70 | *Otherwise, it can also be DifferentialForm::DivConforming, in which case the | ||
| 71 | *normal component is glued continuously. Finally, there is the option of not | ||
| 72 | *gluing the B-splines at all if the parameter is | ||
| 73 | *DifferentialForm::Discontinuous. | ||
| 74 | **/ | ||
| 75 | template <typename Derived> | ||
| 76 | class Glue { | ||
| 77 | std::vector<std::array<int, 4>> edges_; | ||
| 78 | int dofs_after_glue; | ||
| 79 | Eigen::SparseMatrix<double> glue_mat_; | ||
| 80 | |||
| 81 | public: | ||
| 82 | /** | ||
| 83 | * \brief Constructor for the Glue class. | ||
| 84 | * | ||
| 85 | * This constructor initializes a Glue object with the provided SuperSpace and | ||
| 86 | * Projector objects. | ||
| 87 | * | ||
| 88 | * \param superspace The SuperSpace to handle the basis functions. | ||
| 89 | * \param proj The Projector which provides information on the dofs before | ||
| 90 | * gluing. | ||
| 91 | */ | ||
| 92 |
1/2✓ Branch 2 taken 144 times.
✗ Branch 3 not taken.
|
272 | Glue(const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { |
| 93 |
1/2✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
|
272 | init_Glue(superspace, proj); |
| 94 | 272 | } | |
| 95 | /** | ||
| 96 | * \brief Initializes the AnsatzSpace object. | ||
| 97 | * | ||
| 98 | * This function initializes the member variables of the Glue object. | ||
| 99 | * | ||
| 100 | * \param superspace The SuperSpace to handle the basis functions. | ||
| 101 | * \param proj The Projector which provides information on the dofs before | ||
| 102 | * gluing. | ||
| 103 | */ | ||
| 104 | 272 | void init_Glue(const SuperSpace<Derived> &superspace, | |
| 105 | const Projector<Derived> &proj) { | ||
| 106 | 272 | edges_ = superspace.get_mesh().get_element_tree().patchTopologyInfo(); | |
| 107 |
1/2✓ Branch 2 taken 144 times.
✗ Branch 3 not taken.
|
272 | glue_mat_ = assembleGlueMatrix(superspace, proj); |
| 108 | 272 | return; | |
| 109 | } | ||
| 110 | /** | ||
| 111 | * \brief Returns Glue matrix to assemble continuous B-splines. | ||
| 112 | * | ||
| 113 | * This function returns the Glue matrix which transforms smooth B-splines to | ||
| 114 | * continuous shape functions across patch boundaries. If continuity is not | ||
| 115 | * needed this matrix is the identity. Other wise this matrix is a tall thin | ||
| 116 | * matrix combining degrees of freedom. | ||
| 117 | * | ||
| 118 | * \return Eigen::SparseMatrix<double> Glue matrix | ||
| 119 | */ | ||
| 120 | 272 | Eigen::SparseMatrix<double> get_glue_matrix() const { return glue_mat_; } | |
| 121 | /** | ||
| 122 | * \brief Generates a list of degrees of freedom (DOFs) identifications. | ||
| 123 | * | ||
| 124 | * This function creates a list of degrees of freedom (DOFs) identifications | ||
| 125 | * which need to be glued together. The return type is a std::vector of | ||
| 126 | * GlueRoutines::dofIdentification which collects the DOFs to be glued. | ||
| 127 | * Furthermore, the coef in this struct denotes the orientation. | ||
| 128 | * | ||
| 129 | * \param superspace The SuperSpace reference to handle basis functions. | ||
| 130 | * \param proj The Projector contains the number of DOFs. | ||
| 131 | * \return A vector of GlueRoutines::dofIdentification containing the DOF | ||
| 132 | * identifications. | ||
| 133 | */ | ||
| 134 | 272 | inline std::vector<GlueRoutines::dofIdentification> makeDofIdentificationList( | |
| 135 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 136 | return GlueRoutines::glue_identificationmaker_< | ||
| 137 | Derived, | ||
| 138 | 272 | LinearOperatorTraits<Derived>::Form>::makeIdentification(edges_, | |
| 139 | superspace, | ||
| 140 | 272 | proj); | |
| 141 | } | ||
| 142 | /** | ||
| 143 | * \brief Assembles the Glue matrix according to the DOF identification list. | ||
| 144 | * | ||
| 145 | * This function first sorts the DOFs in each identification list. Than the | ||
| 146 | * identification lists get sorted according to the first entries. After that | ||
| 147 | * the DOFs get stored in the Glue matrix. | ||
| 148 | * | ||
| 149 | * \param superspace The SuperSpace reference to handle basis functions. | ||
| 150 | * \param proj The Projector contains the number of DOFs. | ||
| 151 | * \return The Glue matrix. | ||
| 152 | */ | ||
| 153 | 272 | Eigen::SparseMatrix<double> assembleGlueMatrix( | |
| 154 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 155 | 272 | const int pre_dofs = proj.get_dofs_after_projector(); | |
| 156 | // The dofs that need to be identified with each other are divided into | ||
| 157 | // master and slaves, where the master is the firs w.r.t. the tp-ordering | ||
| 158 | 272 | std::vector<Eigen::Triplet<double>> trips; | |
| 159 |
1/2✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
|
272 | std::vector<GlueRoutines::dofIdentification> dof_id = |
| 160 | makeDofIdentificationList(superspace, proj); | ||
| 161 | |||
| 162 | // We keep track of certain information | ||
| 163 |
1/2✓ Branch 2 taken 144 times.
✗ Branch 3 not taken.
|
272 | std::vector<bool> dof_is_slave(pre_dofs, false); |
| 164 |
1/2✓ Branch 2 taken 144 times.
✗ Branch 3 not taken.
|
272 | std::vector<bool> dof_is_master(pre_dofs, false); |
| 165 | |||
| 166 | // Here, we fill the two vectors above | ||
| 167 | 272 | int number_of_slaves = 0; | |
| 168 | 272 | int number_of_boundary_dofs = 0; | |
| 169 |
3/4✓ Branch 4 taken 6812 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 6812 times.
✓ Branch 10 taken 144 times.
|
13804 | for (auto dofset : dof_id) { |
| 170 | // This sorting is required, since by construction only dofs[0]<dofs[1] is | ||
| 171 | // given, but in the case of corners and H1 discretizations, it might | ||
| 172 | // happen that dofs[0]>dofs[2]. Moreover, we count how many dofs we "glue | ||
| 173 | // away" or skipped on the boundary. | ||
| 174 |
1/2✓ Branch 3 taken 6812 times.
✗ Branch 4 not taken.
|
13532 | std::sort(dofset.dofs.begin(), dofset.dofs.end()); |
| 175 |
1/2✓ Branch 2 taken 6812 times.
✗ Branch 3 not taken.
|
13532 | dof_is_master[dofset.dofs[0]] = true; |
| 176 |
2/2✓ Branch 1 taken 4044 times.
✓ Branch 2 taken 2768 times.
|
13532 | if (dofset.dofs.size() == 1) ++number_of_boundary_dofs; |
| 177 |
2/2✓ Branch 1 taken 3040 times.
✓ Branch 2 taken 6812 times.
|
19516 | for (int i = 1; i < dofset.dofs.size(); ++i) { |
| 178 |
1/2✓ Branch 2 taken 3040 times.
✗ Branch 3 not taken.
|
5984 | dof_is_slave[dofset.dofs[i]] = true; |
| 179 | 5984 | ++number_of_slaves; | |
| 180 | } | ||
| 181 | } | ||
| 182 | |||
| 183 | // Now, we sort w.r.t. the masters. | ||
| 184 |
1/2✓ Branch 3 taken 144 times.
✗ Branch 4 not taken.
|
272 | std::sort(dof_id.begin(), dof_id.end(), |
| 185 | 45710 | [](GlueRoutines::dofIdentification a, | |
| 186 | GlueRoutines::dofIdentification b) { | ||
| 187 | 45710 | return a.dofs[0] < b.dofs[0]; | |
| 188 | }); | ||
| 189 | |||
| 190 | 272 | const int post_dofs = pre_dofs - number_of_slaves - number_of_boundary_dofs; | |
| 191 | 272 | const int glued_dofs = pre_dofs - number_of_slaves; | |
| 192 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 144 times.
|
272 | assert(post_dofs != 0 && "All degrees of freedom are on the boundary!"); |
| 193 | // This block is the heart of the algorithms. Skip keeps track of how many | ||
| 194 | // slaves have been skipped already, and master_index keeps track on which | ||
| 195 | // master is about to be glued next. post_index keeps track on how many | ||
| 196 | // dofs are on the boundary of the patches and therefore skipped. | ||
| 197 | 272 | int skip = 0; | |
| 198 | 272 | int master_index = 0; | |
| 199 | 272 | int post_index = 0; | |
| 200 |
2/2✓ Branch 0 taken 43596 times.
✓ Branch 1 taken 144 times.
|
82524 | for (int i = 0; i < glued_dofs; ++i) { |
| 201 |
6/8✓ Branch 1 taken 46602 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3006 times.
✓ Branch 5 taken 43596 times.
✓ Branch 6 taken 3006 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 3006 times.
✓ Branch 9 taken 43596 times.
|
88202 | while (dof_is_slave[i + skip] && ((i + skip) < pre_dofs)) ++skip; |
| 202 | 82252 | const int pre_index = i + skip; | |
| 203 | // skip dofs on patch boundary without a slave patch | ||
| 204 |
7/8✓ Branch 1 taken 43596 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6812 times.
✓ Branch 5 taken 36784 times.
✓ Branch 8 taken 4044 times.
✓ Branch 9 taken 2768 times.
✓ Branch 10 taken 4044 times.
✓ Branch 11 taken 39552 times.
|
82252 | if (dof_is_master[pre_index] && dof_id[master_index].dofs.size() == 1) { |
| 205 | 8076 | ++master_index; | |
| 206 | 8076 | continue; | |
| 207 | } | ||
| 208 |
1/2✓ Branch 2 taken 39552 times.
✗ Branch 3 not taken.
|
74176 | trips.push_back(Eigen::Triplet<double>(pre_index, post_index, 1)); |
| 209 | // The dof cannot be declared slave and master at the same time | ||
| 210 |
5/8✓ Branch 1 taken 39552 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2768 times.
✓ Branch 5 taken 36784 times.
✓ Branch 7 taken 2768 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2768 times.
|
74176 | assert(!(dof_is_master[pre_index] && dof_is_slave[pre_index])); |
| 211 |
3/4✓ Branch 1 taken 39552 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2768 times.
✓ Branch 5 taken 36784 times.
|
74176 | if (dof_is_master[pre_index]) { |
| 212 | // The dofs in dof_id don't know they might be moved forward. So the | ||
| 213 | // smallest dof of those to be identified should coincide with the | ||
| 214 | // pre_index. | ||
| 215 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 2768 times.
|
5456 | assert(pre_index == dof_id[master_index].dofs[0]); |
| 216 | 5456 | const int number_of_partners = dof_id[master_index].dofs.size(); | |
| 217 |
2/2✓ Branch 0 taken 3040 times.
✓ Branch 1 taken 2768 times.
|
11440 | for (int j = 1; j < number_of_partners; ++j) { |
| 218 |
1/2✓ Branch 3 taken 3040 times.
✗ Branch 4 not taken.
|
5984 | trips.push_back(Eigen::Triplet<double>(dof_id[master_index].dofs[j], |
| 219 | post_index, | ||
| 220 | 5984 | dof_id[master_index].coef)); | |
| 221 | } | ||
| 222 | 5456 | ++master_index; | |
| 223 | } | ||
| 224 | 74176 | ++post_index; | |
| 225 | } | ||
| 226 | |||
| 227 |
1/2✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
|
272 | Eigen::SparseMatrix<double> glue_matrix(pre_dofs, post_dofs); |
| 228 |
1/2✓ Branch 3 taken 144 times.
✗ Branch 4 not taken.
|
272 | glue_matrix.setFromTriplets(trips.begin(), trips.end()); |
| 229 | |||
| 230 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 144 times.
|
272 | assert(trips.size() == (pre_dofs - number_of_boundary_dofs)); |
| 231 | 544 | return glue_matrix; | |
| 232 | 272 | } | |
| 233 | }; | ||
| 234 | |||
| 235 | namespace GlueRoutines { | ||
| 236 | |||
| 237 | /** | ||
| 238 | * \brief This routine collects the indices of DOFs with support on an edge. | ||
| 239 | * | ||
| 240 | * The indices are computed depending on the edge case and number of DOFs in the | ||
| 241 | * x and y direction of the tensor product. The shift parameter can be used to | ||
| 242 | * shift all dofs according to the patch index or vector component. | ||
| 243 | * | ||
| 244 | * Edge Case | ||
| 245 | * | ||
| 246 | * 0 = (0,0)->(1,0) | ||
| 247 | * | ||
| 248 | * 1 = (1,0)->(1,1) | ||
| 249 | * | ||
| 250 | * 2 = (1,1)->(0,1) | ||
| 251 | * | ||
| 252 | * 3 = (0,1)->(0,0) | ||
| 253 | * | ||
| 254 | * \param edgeCase Edge case of the patch to collect the DOFs. | ||
| 255 | * \param dimXdir Number of DOFs in the x direction of the tensor product. | ||
| 256 | * \param dimYdir Number of DOFs in the y direction of the tensor product. | ||
| 257 | * \param shift Offset of the DOFs. | ||
| 258 | */ | ||
| 259 | 4224 | std::vector<int> getEdgeDofIndices(int edgeCase, int dimXdir, int dimYdir, | |
| 260 | int shift) { | ||
| 261 | 4224 | std::vector<int> out; | |
| 262 |
5/5✓ Branch 0 taken 669 times.
✓ Branch 1 taken 669 times.
✓ Branch 2 taken 669 times.
✓ Branch 3 taken 669 times.
✓ Branch 4 taken 1548 times.
|
4224 | switch (edgeCase) { |
| 263 | 669 | case (0): { | |
| 264 |
1/2✓ Branch 1 taken 669 times.
✗ Branch 2 not taken.
|
669 | out.reserve(dimXdir); |
| 265 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 669 times.
|
669 | assert( |
| 266 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 267 | dimXdir <= dimYdir); | ||
| 268 |
2/2✓ Branch 0 taken 3739 times.
✓ Branch 1 taken 669 times.
|
4408 | for (int i = 0; i < dimXdir; ++i) { |
| 269 |
1/2✓ Branch 1 taken 3739 times.
✗ Branch 2 not taken.
|
3739 | out.push_back(i + shift); |
| 270 | } | ||
| 271 | 669 | return out; | |
| 272 | }; | ||
| 273 | 669 | case (1): { | |
| 274 |
1/2✓ Branch 1 taken 669 times.
✗ Branch 2 not taken.
|
669 | out.reserve(dimYdir); |
| 275 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 669 times.
|
669 | assert( |
| 276 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 277 | dimYdir <= dimXdir); | ||
| 278 |
2/2✓ Branch 0 taken 3739 times.
✓ Branch 1 taken 669 times.
|
4408 | for (int i = 0; i < dimYdir; ++i) { |
| 279 |
1/2✓ Branch 1 taken 3739 times.
✗ Branch 2 not taken.
|
3739 | out.push_back(dimXdir * (i + 1) - 1 + shift); |
| 280 | } | ||
| 281 | 669 | return out; | |
| 282 | }; | ||
| 283 | 669 | case (2): { | |
| 284 |
1/2✓ Branch 1 taken 669 times.
✗ Branch 2 not taken.
|
669 | out.reserve(dimXdir); |
| 285 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 669 times.
|
669 | assert( |
| 286 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 287 | dimXdir <= dimYdir); | ||
| 288 |
2/2✓ Branch 0 taken 3739 times.
✓ Branch 1 taken 669 times.
|
4408 | for (int i = 0; i < dimXdir; ++i) { |
| 289 |
1/2✓ Branch 1 taken 3739 times.
✗ Branch 2 not taken.
|
3739 | out.push_back(dimXdir * (dimYdir - 1) + i + shift); |
| 290 | } | ||
| 291 | 669 | return out; | |
| 292 | }; | ||
| 293 | 669 | case (3): { | |
| 294 |
1/2✓ Branch 1 taken 669 times.
✗ Branch 2 not taken.
|
669 | out.reserve(dimYdir); |
| 295 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 669 times.
|
669 | assert( |
| 296 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 297 | dimYdir <= dimXdir); | ||
| 298 |
2/2✓ Branch 0 taken 3739 times.
✓ Branch 1 taken 669 times.
|
4408 | for (int i = 0; i < dimYdir; ++i) { |
| 299 |
1/2✓ Branch 1 taken 3739 times.
✗ Branch 2 not taken.
|
3739 | out.push_back(dimXdir * i + shift); |
| 300 | } | ||
| 301 | 669 | return out; | |
| 302 | }; | ||
| 303 | 1548 | default: { | |
| 304 | }; | ||
| 305 | // An edge might have a -1 index. This occurs only when no partner could | ||
| 306 | // be found, and the -1 is the placeholder of the missing partner. | ||
| 307 | } | ||
| 308 | 1548 | return out; | |
| 309 | ✗ | } | |
| 310 | |||
| 311 | // The following 7 functions figure out how edges and the vector components | ||
| 312 | // are oriented w.r.t. each other. | ||
| 313 | 6196 | inline bool edgeIsForwardParametrized(int edgeCase) { | |
| 314 |
4/4✓ Branch 0 taken 5071 times.
✓ Branch 1 taken 1125 times.
✓ Branch 2 taken 1129 times.
✓ Branch 3 taken 3942 times.
|
6196 | return edgeCase == 0 || edgeCase == 1; |
| 315 | } | ||
| 316 | 3028 | inline bool edgeIsBackwardsParametrized(int edgeCase) { | |
| 317 | 3028 | return !(edgeIsForwardParametrized(edgeCase)); | |
| 318 | } | ||
| 319 | 3128 | inline bool normalComponentIsInwardDirected(int edgeCase) { | |
| 320 |
4/4✓ Branch 0 taken 2559 times.
✓ Branch 1 taken 569 times.
✓ Branch 2 taken 571 times.
✓ Branch 3 taken 1988 times.
|
3128 | return edgeCase == 0 || edgeCase == 3; |
| 321 | } | ||
| 322 | 1534 | inline bool normalComponentIsOutwardDirected(int edgeCase) { | |
| 323 | 1534 | return !(normalComponentIsInwardDirected(edgeCase)); | |
| 324 | } | ||
| 325 | |||
| 326 | 2112 | inline bool reverseParametrized(const std::array<int, 4> &edge) { | |
| 327 |
2/2✓ Branch 2 taken 916 times.
✓ Branch 3 taken 140 times.
|
3168 | return ((edgeIsForwardParametrized(edge[2]) && |
| 328 | 1056 | edgeIsForwardParametrized(edge[3])) || | |
| 329 |
2/2✓ Branch 2 taken 914 times.
✓ Branch 3 taken 142 times.
|
3028 | (edgeIsBackwardsParametrized(edge[2]) && |
| 330 | 1056 | edgeIsBackwardsParametrized(edge[3]))) | |
| 331 |
4/4✓ Branch 0 taken 1056 times.
✓ Branch 1 taken 1056 times.
✓ Branch 2 taken 1056 times.
✓ Branch 3 taken 916 times.
|
5280 | ? true |
| 332 | 2112 | : false; | |
| 333 | } | ||
| 334 | |||
| 335 | 1064 | inline int glueCoefficientDivergenceConforming(const std::array<int, 4> &edge) { | |
| 336 |
2/2✓ Branch 2 taken 466 times.
✓ Branch 3 taken 64 times.
|
1594 | return ((normalComponentIsInwardDirected(edge[2]) && |
| 337 |
2/2✓ Branch 2 taken 534 times.
✓ Branch 3 taken 466 times.
|
1530 | normalComponentIsInwardDirected(edge[3])) || |
| 338 |
2/2✓ Branch 2 taken 454 times.
✓ Branch 3 taken 80 times.
|
1534 | (normalComponentIsOutwardDirected(edge[2]) && |
| 339 | 534 | normalComponentIsOutwardDirected(edge[3]))) | |
| 340 |
2/2✓ Branch 0 taken 530 times.
✓ Branch 1 taken 534 times.
|
2128 | ? -1 |
| 341 | 1064 | : 1; | |
| 342 | } | ||
| 343 | 6384 | inline bool edgeToBeGluedInFirstComp(int edgeCase) { | |
| 344 |
4/4✓ Branch 0 taken 5373 times.
✓ Branch 1 taken 1011 times.
✓ Branch 2 taken 4362 times.
✓ Branch 3 taken 1011 times.
|
6384 | return (edgeCase == 0 || edgeCase == 2) ? false : true; |
| 345 | } | ||
| 346 | |||
| 347 | /** | ||
| 348 | * \brief Specification of the struct for the discontinuous case. | ||
| 349 | * | ||
| 350 | * Nothing needs to be glued. | ||
| 351 | */ | ||
| 352 | template <typename Derived> | ||
| 353 | struct glue_identificationmaker_<Derived, DifferentialForm::Discontinuous> { | ||
| 354 | 11 | static std::vector<dofIdentification> makeIdentification( | |
| 355 | const std::vector<std::array<int, 4>> &edges_, | ||
| 356 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 357 | 11 | return {}; | |
| 358 | } | ||
| 359 | }; | ||
| 360 | |||
| 361 | /** | ||
| 362 | * \brief Specification of the struct for the scalar continuous case. | ||
| 363 | * | ||
| 364 | * The scalar case: Here, we have only one vector component to worry about. | ||
| 365 | * However, we need to take care of edges, since here multiple (more than 2) | ||
| 366 | * dofs will be reduced to one. This is also the reason we we work with the | ||
| 367 | * dofIdentification struct and not already with Eigen::Triplet. The function | ||
| 368 | * builds upon the fact, that px == py == p. | ||
| 369 | */ | ||
| 370 | template <typename Derived> | ||
| 371 | struct glue_identificationmaker_<Derived, DifferentialForm::Continuous> { | ||
| 372 | 66 | static std::vector<dofIdentification> makeIdentification( | |
| 373 | const std::vector<std::array<int, 4>> &edges_, | ||
| 374 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 375 | // We check if the space can even be continuous globally. | ||
| 376 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 66 times.
|
66 | assert(superspace.get_polynomial_degree() >= 1); |
| 377 | 66 | const int one_d_dim = superspace.get_polynomial_degree() + 1 + | |
| 378 | 66 | proj.get_knot_repetition() * | |
| 379 | 66 | ((1 << superspace.get_refinement_level()) - 1); | |
| 380 | 66 | const int dofs_per_patch = one_d_dim * one_d_dim; | |
| 381 | |||
| 382 | 66 | std::vector<GlueRoutines::dofIdentification> out; | |
| 383 |
1/2✓ Branch 2 taken 66 times.
✗ Branch 3 not taken.
|
66 | out.reserve(edges_.size() * one_d_dim); |
| 384 | |||
| 385 |
1/2✓ Branch 3 taken 66 times.
✗ Branch 4 not taken.
|
66 | std::vector<int> already_stored_in(proj.get_dofs_after_projector(), -1); |
| 386 | 66 | int d_count = 0; | |
| 387 |
2/2✓ Branch 7 taken 1048 times.
✓ Branch 8 taken 66 times.
|
1114 | for (auto edge : edges_) { |
| 388 | // The shift is to account for different patches, since the | ||
| 389 | // getEdgeDof-routines only can enumerate w.r.t. patch 0. | ||
| 390 | 1048 | const int shift_e1 = edge[0] * dofs_per_patch; | |
| 391 | 1048 | const int shift_e2 = edge[1] * dofs_per_patch; | |
| 392 | 1048 | const bool needReversion = reverseParametrized(edge); | |
| 393 | |||
| 394 |
1/2✓ Branch 1 taken 1048 times.
✗ Branch 2 not taken.
|
1048 | std::vector<int> dofs_e1 = |
| 395 | 1048 | getEdgeDofIndices(edge[2], one_d_dim, one_d_dim, shift_e1); | |
| 396 |
1/2✓ Branch 1 taken 1048 times.
✗ Branch 2 not taken.
|
1048 | std::vector<int> dofs_e2 = |
| 397 | 1048 | getEdgeDofIndices(edge[3], one_d_dim, one_d_dim, shift_e2); | |
| 398 | |||
| 399 | // Check if the edge is hanging. For screens one would need an "else". | ||
| 400 |
5/6✓ Branch 1 taken 280 times.
✓ Branch 2 taken 768 times.
✓ Branch 4 taken 280 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 280 times.
✓ Branch 7 taken 768 times.
|
1048 | if (edge[1] > -1 && edge[0] > -1) { |
| 401 |
2/2✓ Branch 1 taken 1648 times.
✓ Branch 2 taken 280 times.
|
1928 | for (int i = 0; i < one_d_dim; ++i) { |
| 402 | 1648 | GlueRoutines::dofIdentification d; | |
| 403 |
2/2✓ Branch 0 taken 824 times.
✓ Branch 1 taken 824 times.
|
1648 | const int j = needReversion ? one_d_dim - 1 - i : i; |
| 404 | // const int j = i; | ||
| 405 |
1/2✓ Branch 2 taken 1648 times.
✗ Branch 3 not taken.
|
1648 | assert(dofs_e1[i] != dofs_e2[j] && |
| 406 | "If this happens something went horribly wrong."); | ||
| 407 | |||
| 408 | 1648 | const int small_dof = std::min(dofs_e1[i], dofs_e2[j]); | |
| 409 | 1648 | const int large_dof = std::max(dofs_e1[i], dofs_e2[j]); | |
| 410 | |||
| 411 | // In the continuous case, more than two dofs might be glued together. | ||
| 412 | // Therefore, we must check if we know one of the dofs already. If | ||
| 413 | // yes, we just add the new one to the dof_id, if not, we make a new | ||
| 414 | // identification group. | ||
| 415 |
2/2✓ Branch 1 taken 288 times.
✓ Branch 2 taken 1360 times.
|
1648 | if (already_stored_in[small_dof] > -1) { |
| 416 | // It could be that the large_dof is already matched with another | ||
| 417 | // dof, i.e., in the case of a 'circle'. Then there is nothing to | ||
| 418 | // do. If not, we match it with the master of the partner. | ||
| 419 |
2/2✓ Branch 1 taken 272 times.
✓ Branch 2 taken 16 times.
|
288 | if (already_stored_in[large_dof] == -1) { |
| 420 |
1/2✓ Branch 3 taken 272 times.
✗ Branch 4 not taken.
|
272 | out[already_stored_in[small_dof]].dofs.push_back(large_dof); |
| 421 | 272 | already_stored_in[large_dof] = already_stored_in[small_dof]; | |
| 422 | } else { | ||
| 423 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
16 | if (!(already_stored_in[small_dof] == |
| 424 | already_stored_in[large_dof])) { | ||
| 425 | // This case is tricky. Assume that we have to identify four | ||
| 426 | // DOFs with each other, but they have already been assigned in | ||
| 427 | // pairs of two. Then we need to reverse this process. First, we | ||
| 428 | // grab the two storage locations. | ||
| 429 | ✗ | const int small_store = std::min(already_stored_in[large_dof], | |
| 430 | already_stored_in[small_dof]); | ||
| 431 | ✗ | const int large_store = std::max(already_stored_in[large_dof], | |
| 432 | already_stored_in[small_dof]); | ||
| 433 | ✗ | for (auto dfs : out[large_store].dofs) { | |
| 434 | // now we put all of those in the larger location into the | ||
| 435 | // smaller one. | ||
| 436 | ✗ | already_stored_in[dfs] = small_store; | |
| 437 | ✗ | for (auto other_dfs : out[small_store].dofs) { | |
| 438 | ✗ | assert(dfs != other_dfs); | |
| 439 | } | ||
| 440 | ✗ | out[small_store].dofs.push_back(dfs); | |
| 441 | } | ||
| 442 | // Now we set the larger storage location to empty. | ||
| 443 | ✗ | out[large_store].dofs = {}; | |
| 444 | ✗ | out[large_store].coef = 0; | |
| 445 | } | ||
| 446 | } | ||
| 447 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1360 times.
|
1360 | } else if (already_stored_in[large_dof] > -1) { |
| 448 | // It could be that the small_dof is already matched with another | ||
| 449 | // dof, i.e., in the case of a 'circle'. Then there is nothing to | ||
| 450 | // do. If not, we match it with the master of the partner. | ||
| 451 | ✗ | if (already_stored_in[small_dof] == -1) { | |
| 452 | ✗ | out[already_stored_in[large_dof]].dofs.push_back(small_dof); | |
| 453 | ✗ | already_stored_in[small_dof] = already_stored_in[large_dof]; | |
| 454 | } else { | ||
| 455 | ✗ | if (!(already_stored_in[small_dof] == | |
| 456 | already_stored_in[large_dof])) { | ||
| 457 | // This case is tricky. Assume that we have to identify four | ||
| 458 | // DOFs with each other, but they have already been assigned in | ||
| 459 | // pairs of two. Then we need to reverse this process. First, we | ||
| 460 | // grab the two storage locations. | ||
| 461 | ✗ | const int small_store = std::min(already_stored_in[large_dof], | |
| 462 | already_stored_in[small_dof]); | ||
| 463 | ✗ | const int large_store = std::max(already_stored_in[large_dof], | |
| 464 | already_stored_in[small_dof]); | ||
| 465 | ✗ | for (auto dfs : out[large_store].dofs) { | |
| 466 | // now we put all of those in the larger location into the | ||
| 467 | // smaller one. | ||
| 468 | ✗ | already_stored_in[dfs] = small_store; | |
| 469 | ✗ | for (auto other_dfs : out[small_store].dofs) { | |
| 470 | ✗ | assert(dfs != other_dfs); | |
| 471 | } | ||
| 472 | ✗ | out[small_store].dofs.push_back(dfs); | |
| 473 | } | ||
| 474 | // Now we set the larger storage location to empty. | ||
| 475 | ✗ | out[large_store].dofs = {}; | |
| 476 | ✗ | out[large_store].coef = 0; | |
| 477 | } | ||
| 478 | } | ||
| 479 | } else { | ||
| 480 | // With the exception of corners, this will be the default case. | ||
| 481 | // We just add a pair of dofs to the bookkeeping. | ||
| 482 |
1/2✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
|
1360 | d.dofs.push_back(small_dof); |
| 483 |
1/2✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
|
1360 | d.dofs.push_back(large_dof); |
| 484 | 1360 | already_stored_in[small_dof] = d_count; | |
| 485 | 1360 | already_stored_in[large_dof] = d_count; | |
| 486 | 1360 | d_count++; | |
| 487 | 1360 | d.coef = 1; | |
| 488 |
1/2✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
|
1360 | out.push_back(d); |
| 489 | } | ||
| 490 | } | ||
| 491 | } | ||
| 492 | } | ||
| 493 | // Now we need to clean up the empty dofsets, since subsequent routines | ||
| 494 | // assume there to be at least one element. | ||
| 495 |
2/2✓ Branch 4 taken 1360 times.
✓ Branch 5 taken 66 times.
|
1426 | for (auto x = out.begin(); x != out.end(); ++x) { |
| 496 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 1360 times.
|
1360 | if ((*x).dofs.size() == 0) { |
| 497 | ✗ | out.erase(x); | |
| 498 | } | ||
| 499 | } | ||
| 500 | 132 | return out; | |
| 501 | 66 | } | |
| 502 | }; | ||
| 503 | |||
| 504 | /** | ||
| 505 | * \brief Specification of the struct for the vector Div conforming case. | ||
| 506 | * | ||
| 507 | * The Maxwell case: Here, the identification is 1-to-1, but we need to take | ||
| 508 | * care to glue the right vector component. The function builds upon the fact, | ||
| 509 | * that px == py == p. | ||
| 510 | */ | ||
| 511 | template <typename Derived> | ||
| 512 | struct glue_identificationmaker_<Derived, DifferentialForm::DivConforming> { | ||
| 513 | 67 | static std::vector<dofIdentification> makeIdentification( | |
| 514 | const std::vector<std::array<int, 4>> &edges_, | ||
| 515 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 516 | // since we assume px == py and uniform refinement in both directions, there | ||
| 517 | // will be a small_dim and a large_dim in every vector component. | ||
| 518 | 67 | const int small_dim = superspace.get_polynomial_degree() + | |
| 519 | 67 | proj.get_knot_repetition() * | |
| 520 | 67 | ((1 << superspace.get_refinement_level()) - 1); | |
| 521 | 67 | const int large_dim = superspace.get_polynomial_degree() + 1 + | |
| 522 | 67 | proj.get_knot_repetition() * | |
| 523 | 67 | ((1 << superspace.get_refinement_level()) - 1); | |
| 524 | 67 | const int dofs_per_patch_per_component = small_dim * large_dim; | |
| 525 | // sanity check | ||
| 526 |
1/2✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
|
67 | assert(dofs_per_patch_per_component == |
| 527 | (proj.get_dofs_after_projector() / | ||
| 528 | (superspace.get_number_of_patches() * 2)) && | ||
| 529 | "The assembly of the glue matrix is highly specific to the space. " | ||
| 530 | "Something went wrong; is the discrete space correct?"); | ||
| 531 | |||
| 532 | 67 | std::vector<GlueRoutines::dofIdentification> out; | |
| 533 |
1/2✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
|
67 | out.reserve(edges_.size() * small_dim); |
| 534 | |||
| 535 |
2/2✓ Branch 7 taken 1064 times.
✓ Branch 8 taken 67 times.
|
1131 | for (auto edge : edges_) { |
| 536 |
3/4✓ Branch 2 taken 780 times.
✓ Branch 3 taken 284 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 780 times.
|
1064 | assert(edge[0] <= edge[1] || edge[1] == -1); |
| 537 | // The DOFs that need to be glued on edges 0 and 2 are in the second | ||
| 538 | // vector component. The remainder of the dof-index-shift is to account | ||
| 539 | // for dofs on other patches. | ||
| 540 | 1064 | const int shift_e1 = edge[0] * dofs_per_patch_per_component + | |
| 541 | 1064 | (edgeToBeGluedInFirstComp(edge[2]) | |
| 542 |
2/2✓ Branch 0 taken 534 times.
✓ Branch 1 taken 530 times.
|
1594 | ? 0 |
| 543 | : (dofs_per_patch_per_component * | ||
| 544 | 530 | superspace.get_number_of_patches())); | |
| 545 | 1064 | const int shift_e2 = edge[1] * dofs_per_patch_per_component + | |
| 546 | 1064 | (edgeToBeGluedInFirstComp(edge[3]) | |
| 547 |
2/2✓ Branch 0 taken 920 times.
✓ Branch 1 taken 144 times.
|
1208 | ? 0 |
| 548 | : (dofs_per_patch_per_component * | ||
| 549 | 144 | superspace.get_number_of_patches())); | |
| 550 | |||
| 551 | 1064 | const int xdir_dim_1 = | |
| 552 |
2/2✓ Branch 2 taken 534 times.
✓ Branch 3 taken 530 times.
|
1064 | edgeToBeGluedInFirstComp(edge[2]) ? large_dim : small_dim; |
| 553 | 1064 | const int ydir_dim_1 = | |
| 554 |
2/2✓ Branch 2 taken 534 times.
✓ Branch 3 taken 530 times.
|
1064 | edgeToBeGluedInFirstComp(edge[2]) ? small_dim : large_dim; |
| 555 | 1064 | const int xdir_dim_2 = | |
| 556 |
2/2✓ Branch 2 taken 920 times.
✓ Branch 3 taken 144 times.
|
1064 | edgeToBeGluedInFirstComp(edge[3]) ? large_dim : small_dim; |
| 557 | 1064 | const int ydir_dim_2 = | |
| 558 |
2/2✓ Branch 2 taken 920 times.
✓ Branch 3 taken 144 times.
|
1064 | edgeToBeGluedInFirstComp(edge[3]) ? small_dim : large_dim; |
| 559 |
1/2✓ Branch 1 taken 1064 times.
✗ Branch 2 not taken.
|
1064 | std::vector<int> dofs_e1 = |
| 560 | 1064 | getEdgeDofIndices(edge[2], xdir_dim_1, ydir_dim_1, shift_e1); | |
| 561 |
1/2✓ Branch 1 taken 1064 times.
✗ Branch 2 not taken.
|
1064 | std::vector<int> dofs_e2 = |
| 562 | 1064 | getEdgeDofIndices(edge[3], xdir_dim_2, ydir_dim_2, shift_e2); | |
| 563 | |||
| 564 | 1064 | const int size_of_edge_dofs = dofs_e1.size(); | |
| 565 |
3/4✓ Branch 1 taken 780 times.
✓ Branch 2 taken 284 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 780 times.
|
1064 | assert(size_of_edge_dofs == dofs_e2.size() || edge[1] == -1); |
| 566 | |||
| 567 | 1064 | const bool needReversion = reverseParametrized(edge); | |
| 568 | 1064 | const int coef = glueCoefficientDivergenceConforming(edge); | |
| 569 | |||
| 570 | // Check if the edge is hanging. | ||
| 571 |
5/6✓ Branch 1 taken 284 times.
✓ Branch 2 taken 780 times.
✓ Branch 4 taken 284 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 284 times.
✓ Branch 7 taken 780 times.
|
1064 | if (edge[1] > -1 && edge[0] > -1) { |
| 572 |
2/2✓ Branch 1 taken 1408 times.
✓ Branch 2 taken 284 times.
|
1692 | for (int i = 0; i < size_of_edge_dofs; ++i) { |
| 573 | 1408 | GlueRoutines::dofIdentification d; | |
| 574 |
2/2✓ Branch 0 taken 702 times.
✓ Branch 1 taken 706 times.
|
1408 | const int j = needReversion ? size_of_edge_dofs - 1 - i : i; |
| 575 |
1/2✓ Branch 2 taken 1408 times.
✗ Branch 3 not taken.
|
1408 | assert(dofs_e1[i] != dofs_e2[j] && |
| 576 | "If this happens something went horribly wrong."); | ||
| 577 |
1/2✓ Branch 4 taken 1408 times.
✗ Branch 5 not taken.
|
1408 | d.dofs.push_back(std::min(dofs_e1[i], dofs_e2[j])); |
| 578 |
1/2✓ Branch 4 taken 1408 times.
✗ Branch 5 not taken.
|
1408 | d.dofs.push_back(std::max(dofs_e1[i], dofs_e2[j])); |
| 579 | 1408 | d.coef = coef; | |
| 580 |
1/2✓ Branch 1 taken 1408 times.
✗ Branch 2 not taken.
|
1408 | out.push_back(d); |
| 581 | } | ||
| 582 | } else { | ||
| 583 |
2/2✓ Branch 1 taken 4044 times.
✓ Branch 2 taken 780 times.
|
4824 | for (int i = 0; i < size_of_edge_dofs; ++i) { |
| 584 | // Add only the master dof which is skipped in subsequent routines | ||
| 585 | 4044 | GlueRoutines::dofIdentification d; | |
| 586 |
1/2✓ Branch 2 taken 4044 times.
✗ Branch 3 not taken.
|
4044 | d.dofs.push_back(dofs_e1[i]); |
| 587 |
1/2✓ Branch 1 taken 4044 times.
✗ Branch 2 not taken.
|
4044 | out.push_back(d); |
| 588 | } | ||
| 589 | } | ||
| 590 | } | ||
| 591 |
1/2✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
|
67 | out.shrink_to_fit(); |
| 592 | 67 | return out; | |
| 593 | ✗ | } | |
| 594 | }; | ||
| 595 | }; // namespace GlueRoutines | ||
| 596 | |||
| 597 | } // namespace Bembel | ||
| 598 | #endif // BEMBEL_SRC_ANSATZSPACE_GLUE_HPP_ | ||
| 599 |