| 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 |
3/4✓ Branch 4 taken 2764 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 2764 times.
✓ Branch 10 taken 144 times.
|
5724 | for (auto dofset : dof_id) { |
| 169 | // This sorting is required, since by construction only dofs[0]<dofs[1] is | ||
| 170 | // given, but in the case of corners and H1 discretizations, it might | ||
| 171 | // happen that dofs[0]>dofs[2]. Moreover, we count how many dofs we "glue | ||
| 172 | // away". | ||
| 173 |
1/2✓ Branch 3 taken 2764 times.
✗ Branch 4 not taken.
|
5452 | std::sort(dofset.dofs.begin(), dofset.dofs.end()); |
| 174 |
1/2✓ Branch 2 taken 2764 times.
✗ Branch 3 not taken.
|
5452 | dof_is_master[dofset.dofs[0]] = true; |
| 175 |
2/2✓ Branch 1 taken 3036 times.
✓ Branch 2 taken 2764 times.
|
11432 | for (int i = 1; i < dofset.dofs.size(); ++i) { |
| 176 |
1/2✓ Branch 2 taken 3036 times.
✗ Branch 3 not taken.
|
5980 | dof_is_slave[dofset.dofs[i]] = true; |
| 177 | 5980 | ++number_of_slaves; | |
| 178 | } | ||
| 179 | } | ||
| 180 | |||
| 181 | // Now, we sort w.r.t. the masters. | ||
| 182 |
1/2✓ Branch 3 taken 144 times.
✗ Branch 4 not taken.
|
272 | std::sort(dof_id.begin(), dof_id.end(), |
| 183 | 14300 | [](GlueRoutines::dofIdentification a, | |
| 184 | GlueRoutines::dofIdentification b) { | ||
| 185 | 14300 | return a.dofs[0] < b.dofs[0]; | |
| 186 | }); | ||
| 187 | |||
| 188 | 272 | const int post_dofs = pre_dofs - number_of_slaves; | |
| 189 | |||
| 190 | // This block is the heart of the algorithms. Skip keeps track of how many | ||
| 191 | // slaves have been skipped already, and master_index keeps track on which | ||
| 192 | // master is about to be glued next. | ||
| 193 | 272 | int skip = 0; | |
| 194 | 272 | int master_index = 0; | |
| 195 |
2/2✓ Branch 0 taken 43584 times.
✓ Branch 1 taken 144 times.
|
82512 | for (int i = 0; i < post_dofs; ++i) { |
| 196 | 82240 | const int post_index = i; | |
| 197 |
6/8✓ Branch 1 taken 46586 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3002 times.
✓ Branch 5 taken 43584 times.
✓ Branch 6 taken 3002 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 3002 times.
✓ Branch 9 taken 43584 times.
|
88186 | while (dof_is_slave[i + skip] && ((i + skip) < pre_dofs)) ++skip; |
| 198 | 82240 | const int pre_index = i + skip; | |
| 199 |
1/2✓ Branch 2 taken 43584 times.
✗ Branch 3 not taken.
|
82240 | trips.push_back(Eigen::Triplet<double>(pre_index, post_index, 1)); |
| 200 | // The dof cannot be declared slave and master at the same time | ||
| 201 |
5/8✓ Branch 1 taken 43584 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2764 times.
✓ Branch 5 taken 40820 times.
✓ Branch 7 taken 2764 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2764 times.
|
82240 | assert(!(dof_is_master[pre_index] && dof_is_slave[pre_index])); |
| 202 |
3/4✓ Branch 1 taken 43584 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2764 times.
✓ Branch 5 taken 40820 times.
|
82240 | if (dof_is_master[pre_index]) { |
| 203 | // The dofs in dof_id don't know they might be moved forward. So the | ||
| 204 | // smallest dof of those to be identified should coincide with the | ||
| 205 | // pre_index. | ||
| 206 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 2764 times.
|
5452 | assert(pre_index == dof_id[master_index].dofs[0]); |
| 207 | 5452 | const int number_of_partners = dof_id[master_index].dofs.size(); | |
| 208 |
2/2✓ Branch 0 taken 3036 times.
✓ Branch 1 taken 2764 times.
|
11432 | for (int j = 1; j < number_of_partners; ++j) { |
| 209 |
1/2✓ Branch 3 taken 3036 times.
✗ Branch 4 not taken.
|
5980 | trips.push_back(Eigen::Triplet<double>(dof_id[master_index].dofs[j], |
| 210 | post_index, | ||
| 211 | 5980 | dof_id[master_index].coef)); | |
| 212 | } | ||
| 213 | 5452 | ++master_index; | |
| 214 | } | ||
| 215 | } | ||
| 216 | |||
| 217 |
1/2✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
|
272 | Eigen::SparseMatrix<double> glue_matrix(pre_dofs, post_dofs); |
| 218 |
1/2✓ Branch 3 taken 144 times.
✗ Branch 4 not taken.
|
272 | glue_matrix.setFromTriplets(trips.begin(), trips.end()); |
| 219 | |||
| 220 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 144 times.
|
272 | assert(trips.size() == pre_dofs); |
| 221 | 544 | return glue_matrix; | |
| 222 | 272 | } | |
| 223 | }; | ||
| 224 | |||
| 225 | namespace GlueRoutines { | ||
| 226 | |||
| 227 | /** | ||
| 228 | * \brief This routine collects the indices of DOFs with support on an edge. | ||
| 229 | * | ||
| 230 | * The indices are computed depending on the edge case and number of DOFs in the | ||
| 231 | * x and y direction of the tensor product. The shift parameter can be used to | ||
| 232 | * shift all dofs according to the patch index or vector component. | ||
| 233 | * | ||
| 234 | * Edge Case | ||
| 235 | * | ||
| 236 | * 0 = (0,0)->(1,0) | ||
| 237 | * | ||
| 238 | * 1 = (1,0)->(1,1) | ||
| 239 | * | ||
| 240 | * 2 = (1,1)->(0,1) | ||
| 241 | * | ||
| 242 | * 3 = (0,1)->(0,0) | ||
| 243 | * | ||
| 244 | * \param edgeCase Edge case of the patch to collect the DOFs. | ||
| 245 | * \param dimXdir Number of DOFs in the x direction of the tensor product. | ||
| 246 | * \param dimYdir Number of DOFs in the y direction of the tensor product. | ||
| 247 | * \param shift Offset of the DOFs. | ||
| 248 | */ | ||
| 249 | 4200 | std::vector<int> getEdgeDofIndices(int edgeCase, int dimXdir, int dimYdir, | |
| 250 | int shift) { | ||
| 251 | 4200 | std::vector<int> out; | |
| 252 |
5/5✓ Branch 0 taken 665 times.
✓ Branch 1 taken 665 times.
✓ Branch 2 taken 665 times.
✓ Branch 3 taken 665 times.
✓ Branch 4 taken 1540 times.
|
4200 | switch (edgeCase) { |
| 253 | 665 | case (0): { | |
| 254 |
1/2✓ Branch 1 taken 665 times.
✗ Branch 2 not taken.
|
665 | out.reserve(dimXdir); |
| 255 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 665 times.
|
665 | assert( |
| 256 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 257 | dimXdir <= dimYdir); | ||
| 258 |
2/2✓ Branch 0 taken 3735 times.
✓ Branch 1 taken 665 times.
|
4400 | for (int i = 0; i < dimXdir; ++i) { |
| 259 |
1/2✓ Branch 1 taken 3735 times.
✗ Branch 2 not taken.
|
3735 | out.push_back(i + shift); |
| 260 | } | ||
| 261 | 665 | return out; | |
| 262 | }; | ||
| 263 | 665 | case (1): { | |
| 264 |
1/2✓ Branch 1 taken 665 times.
✗ Branch 2 not taken.
|
665 | out.reserve(dimYdir); |
| 265 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 665 times.
|
665 | assert( |
| 266 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 267 | dimYdir <= dimXdir); | ||
| 268 |
2/2✓ Branch 0 taken 3735 times.
✓ Branch 1 taken 665 times.
|
4400 | for (int i = 0; i < dimYdir; ++i) { |
| 269 |
1/2✓ Branch 1 taken 3735 times.
✗ Branch 2 not taken.
|
3735 | out.push_back(dimXdir * (i + 1) - 1 + shift); |
| 270 | } | ||
| 271 | 665 | return out; | |
| 272 | }; | ||
| 273 | 665 | case (2): { | |
| 274 |
1/2✓ Branch 1 taken 665 times.
✗ Branch 2 not taken.
|
665 | out.reserve(dimXdir); |
| 275 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 665 times.
|
665 | assert( |
| 276 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 277 | dimXdir <= dimYdir); | ||
| 278 |
2/2✓ Branch 0 taken 3735 times.
✓ Branch 1 taken 665 times.
|
4400 | for (int i = 0; i < dimXdir; ++i) { |
| 279 |
1/2✓ Branch 1 taken 3735 times.
✗ Branch 2 not taken.
|
3735 | out.push_back(dimXdir * (dimYdir - 1) + i + shift); |
| 280 | } | ||
| 281 | 665 | return out; | |
| 282 | }; | ||
| 283 | 665 | case (3): { | |
| 284 |
1/2✓ Branch 1 taken 665 times.
✗ Branch 2 not taken.
|
665 | out.reserve(dimYdir); |
| 285 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 665 times.
|
665 | assert( |
| 286 | "This should be a given. If not something went wrong in Glue.hpp" && | ||
| 287 | dimYdir <= dimXdir); | ||
| 288 |
2/2✓ Branch 0 taken 3735 times.
✓ Branch 1 taken 665 times.
|
4400 | for (int i = 0; i < dimYdir; ++i) { |
| 289 |
1/2✓ Branch 1 taken 3735 times.
✗ Branch 2 not taken.
|
3735 | out.push_back(dimXdir * i + shift); |
| 290 | } | ||
| 291 | 665 | return out; | |
| 292 | }; | ||
| 293 | 1540 | default: { | |
| 294 | }; | ||
| 295 | // An edge might have a -1 index. This occurs only when no partner could | ||
| 296 | // be found, and the -1 is the placeholder of the missing partner. | ||
| 297 | } | ||
| 298 | 1540 | return out; | |
| 299 | ✗ | } | |
| 300 | |||
| 301 | // The following 7 functions figure out how edges and the vector components | ||
| 302 | // are oriented w.r.t. each other. | ||
| 303 | 6160 | inline bool edgeIsForwardParametrized(int edgeCase) { | |
| 304 |
4/4✓ Branch 0 taken 5042 times.
✓ Branch 1 taken 1118 times.
✓ Branch 2 taken 1122 times.
✓ Branch 3 taken 3920 times.
|
6160 | return edgeCase == 0 || edgeCase == 1; |
| 305 | } | ||
| 306 | 3010 | inline bool edgeIsBackwardsParametrized(int edgeCase) { | |
| 307 | 3010 | return !(edgeIsForwardParametrized(edgeCase)); | |
| 308 | } | ||
| 309 | 3092 | inline bool normalComponentIsInwardDirected(int edgeCase) { | |
| 310 |
4/4✓ Branch 0 taken 2530 times.
✓ Branch 1 taken 562 times.
✓ Branch 2 taken 564 times.
✓ Branch 3 taken 1966 times.
|
3092 | return edgeCase == 0 || edgeCase == 3; |
| 311 | } | ||
| 312 | 1516 | inline bool normalComponentIsOutwardDirected(int edgeCase) { | |
| 313 | 1516 | return !(normalComponentIsInwardDirected(edgeCase)); | |
| 314 | } | ||
| 315 | |||
| 316 | 2100 | inline bool reverseParametrized(const std::array<int, 4> &edge) { | |
| 317 |
2/2✓ Branch 2 taken 910 times.
✓ Branch 3 taken 140 times.
|
3150 | return ((edgeIsForwardParametrized(edge[2]) && |
| 318 | 1050 | edgeIsForwardParametrized(edge[3])) || | |
| 319 |
2/2✓ Branch 2 taken 910 times.
✓ Branch 3 taken 140 times.
|
3010 | (edgeIsBackwardsParametrized(edge[2]) && |
| 320 | 1050 | edgeIsBackwardsParametrized(edge[3]))) | |
| 321 |
4/4✓ Branch 0 taken 1050 times.
✓ Branch 1 taken 1050 times.
✓ Branch 2 taken 1050 times.
✓ Branch 3 taken 910 times.
|
5250 | ? true |
| 322 | 2100 | : false; | |
| 323 | } | ||
| 324 | |||
| 325 | 1052 | inline int glueCoefficientDivergenceConforming(const std::array<int, 4> &edge) { | |
| 326 |
2/2✓ Branch 2 taken 460 times.
✓ Branch 3 taken 64 times.
|
1576 | return ((normalComponentIsInwardDirected(edge[2]) && |
| 327 |
2/2✓ Branch 2 taken 528 times.
✓ Branch 3 taken 460 times.
|
1512 | normalComponentIsInwardDirected(edge[3])) || |
| 328 |
2/2✓ Branch 2 taken 450 times.
✓ Branch 3 taken 78 times.
|
1516 | (normalComponentIsOutwardDirected(edge[2]) && |
| 329 | 528 | normalComponentIsOutwardDirected(edge[3]))) | |
| 330 |
2/2✓ Branch 0 taken 524 times.
✓ Branch 1 taken 528 times.
|
2104 | ? -1 |
| 331 | 1052 | : 1; | |
| 332 | } | ||
| 333 | 6312 | inline bool edgeToBeGluedInFirstComp(int edgeCase) { | |
| 334 |
4/4✓ Branch 0 taken 5313 times.
✓ Branch 1 taken 999 times.
✓ Branch 2 taken 4314 times.
✓ Branch 3 taken 999 times.
|
6312 | return (edgeCase == 0 || edgeCase == 2) ? false : true; |
| 335 | } | ||
| 336 | |||
| 337 | /** | ||
| 338 | * \brief Specification of the struct for the discontinuous case. | ||
| 339 | * | ||
| 340 | * Nothing needs to be glued. | ||
| 341 | */ | ||
| 342 | template <typename Derived> | ||
| 343 | struct glue_identificationmaker_<Derived, DifferentialForm::Discontinuous> { | ||
| 344 | 11 | static std::vector<dofIdentification> makeIdentification( | |
| 345 | const std::vector<std::array<int, 4>> &edges_, | ||
| 346 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 347 | 11 | return {}; | |
| 348 | } | ||
| 349 | }; | ||
| 350 | |||
| 351 | /** | ||
| 352 | * \brief Specification of the struct for the scalar continuous case. | ||
| 353 | * | ||
| 354 | * The scalar case: Here, we have only one vector component to worry about. | ||
| 355 | * However, we need to take care of edges, since here multiple (more than 2) | ||
| 356 | * dofs will be reduced to one. This is also the reason we we work with the | ||
| 357 | * dofIdentification struct and not already with Eigen::Triplet. The function | ||
| 358 | * builds upon the fact, that px == py == p. | ||
| 359 | */ | ||
| 360 | template <typename Derived> | ||
| 361 | struct glue_identificationmaker_<Derived, DifferentialForm::Continuous> { | ||
| 362 | 66 | static std::vector<dofIdentification> makeIdentification( | |
| 363 | const std::vector<std::array<int, 4>> &edges_, | ||
| 364 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 365 | // We check if the space can even be continuous globally. | ||
| 366 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 66 times.
|
66 | assert(superspace.get_polynomial_degree() >= 1); |
| 367 | 66 | const int one_d_dim = superspace.get_polynomial_degree() + 1 + | |
| 368 | 66 | proj.get_knot_repetition() * | |
| 369 | 66 | ((1 << superspace.get_refinement_level()) - 1); | |
| 370 | 66 | const int dofs_per_patch = one_d_dim * one_d_dim; | |
| 371 | |||
| 372 | 66 | std::vector<GlueRoutines::dofIdentification> out; | |
| 373 |
1/2✓ Branch 2 taken 66 times.
✗ Branch 3 not taken.
|
66 | out.reserve(edges_.size() * one_d_dim); |
| 374 | |||
| 375 |
1/2✓ Branch 3 taken 66 times.
✗ Branch 4 not taken.
|
66 | std::vector<int> already_stored_in(proj.get_dofs_after_projector(), -1); |
| 376 | 66 | int d_count = 0; | |
| 377 |
2/2✓ Branch 7 taken 1048 times.
✓ Branch 8 taken 66 times.
|
1114 | for (auto edge : edges_) { |
| 378 | // The shift is to account for different patches, since the | ||
| 379 | // getEdgeDof-routines only can enumerate w.r.t. patch 0. | ||
| 380 | 1048 | const int shift_e1 = edge[0] * dofs_per_patch; | |
| 381 | 1048 | const int shift_e2 = edge[1] * dofs_per_patch; | |
| 382 | 1048 | const bool needReversion = reverseParametrized(edge); | |
| 383 | |||
| 384 |
1/2✓ Branch 1 taken 1048 times.
✗ Branch 2 not taken.
|
1048 | std::vector<int> dofs_e1 = |
| 385 | 1048 | getEdgeDofIndices(edge[2], one_d_dim, one_d_dim, shift_e1); | |
| 386 |
1/2✓ Branch 1 taken 1048 times.
✗ Branch 2 not taken.
|
1048 | std::vector<int> dofs_e2 = |
| 387 | 1048 | getEdgeDofIndices(edge[3], one_d_dim, one_d_dim, shift_e2); | |
| 388 | |||
| 389 | // Check if the edge is hanging. For screens one would need an "else". | ||
| 390 |
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) { |
| 391 |
2/2✓ Branch 1 taken 1648 times.
✓ Branch 2 taken 280 times.
|
1928 | for (int i = 0; i < one_d_dim; ++i) { |
| 392 | 1648 | GlueRoutines::dofIdentification d; | |
| 393 |
2/2✓ Branch 0 taken 824 times.
✓ Branch 1 taken 824 times.
|
1648 | const int j = needReversion ? one_d_dim - 1 - i : i; |
| 394 | // const int j = i; | ||
| 395 |
1/2✓ Branch 2 taken 1648 times.
✗ Branch 3 not taken.
|
1648 | assert(dofs_e1[i] != dofs_e2[j] && |
| 396 | "If this happens something went horribly wrong."); | ||
| 397 | |||
| 398 | 1648 | const int small_dof = std::min(dofs_e1[i], dofs_e2[j]); | |
| 399 | 1648 | const int large_dof = std::max(dofs_e1[i], dofs_e2[j]); | |
| 400 | |||
| 401 | // In the continuous case, more than two dofs might be glued together. | ||
| 402 | // Therefore, we must check if we know one of the dofs already. If | ||
| 403 | // yes, we just add the new one to the dof_id, if not, we make a new | ||
| 404 | // identification group. | ||
| 405 |
2/2✓ Branch 1 taken 288 times.
✓ Branch 2 taken 1360 times.
|
1648 | if (already_stored_in[small_dof] > -1) { |
| 406 | // It could be that the large_dof is already matched with another | ||
| 407 | // dof, i.e., in the case of a 'circle'. Then there is nothing to | ||
| 408 | // do. If not, we match it with the master of the partner. | ||
| 409 |
2/2✓ Branch 1 taken 272 times.
✓ Branch 2 taken 16 times.
|
288 | if (already_stored_in[large_dof] == -1) { |
| 410 |
1/2✓ Branch 3 taken 272 times.
✗ Branch 4 not taken.
|
272 | out[already_stored_in[small_dof]].dofs.push_back(large_dof); |
| 411 | 272 | already_stored_in[large_dof] = already_stored_in[small_dof]; | |
| 412 | } else { | ||
| 413 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
16 | if (!(already_stored_in[small_dof] == |
| 414 | already_stored_in[large_dof])) { | ||
| 415 | // This case is tricky. Assume that we have to identify four | ||
| 416 | // DOFs with each other, but they have already been assigned in | ||
| 417 | // pairs of two. Then we need to reverse this process. First, we | ||
| 418 | // grab the two storage locations. | ||
| 419 | ✗ | const int small_store = std::min(already_stored_in[large_dof], | |
| 420 | already_stored_in[small_dof]); | ||
| 421 | ✗ | const int large_store = std::max(already_stored_in[large_dof], | |
| 422 | already_stored_in[small_dof]); | ||
| 423 | ✗ | for (auto dfs : out[large_store].dofs) { | |
| 424 | // now we put all of those in the larger location into the | ||
| 425 | // smaller one. | ||
| 426 | ✗ | already_stored_in[dfs] = small_store; | |
| 427 | ✗ | for (auto other_dfs : out[small_store].dofs) { | |
| 428 | ✗ | assert(dfs != other_dfs); | |
| 429 | } | ||
| 430 | ✗ | out[small_store].dofs.push_back(dfs); | |
| 431 | } | ||
| 432 | // Now we set the larger storage location to empty. | ||
| 433 | ✗ | out[large_store].dofs = {}; | |
| 434 | ✗ | out[large_store].coef = 0; | |
| 435 | } | ||
| 436 | } | ||
| 437 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1360 times.
|
1360 | } else if (already_stored_in[large_dof] > -1) { |
| 438 | // It could be that the small_dof is already matched with another | ||
| 439 | // dof, i.e., in the case of a 'circle'. Then there is nothing to | ||
| 440 | // do. If not, we match it with the master of the partner. | ||
| 441 | ✗ | if (already_stored_in[small_dof] == -1) { | |
| 442 | ✗ | out[already_stored_in[large_dof]].dofs.push_back(small_dof); | |
| 443 | ✗ | already_stored_in[small_dof] = already_stored_in[large_dof]; | |
| 444 | } else { | ||
| 445 | ✗ | if (!(already_stored_in[small_dof] == | |
| 446 | already_stored_in[large_dof])) { | ||
| 447 | // This case is tricky. Assume that we have to identify four | ||
| 448 | // DOFs with each other, but they have already been assigned in | ||
| 449 | // pairs of two. Then we need to reverse this process. First, we | ||
| 450 | // grab the two storage locations. | ||
| 451 | ✗ | const int small_store = std::min(already_stored_in[large_dof], | |
| 452 | already_stored_in[small_dof]); | ||
| 453 | ✗ | const int large_store = std::max(already_stored_in[large_dof], | |
| 454 | already_stored_in[small_dof]); | ||
| 455 | ✗ | for (auto dfs : out[large_store].dofs) { | |
| 456 | // now we put all of those in the larger location into the | ||
| 457 | // smaller one. | ||
| 458 | ✗ | already_stored_in[dfs] = small_store; | |
| 459 | ✗ | for (auto other_dfs : out[small_store].dofs) { | |
| 460 | ✗ | assert(dfs != other_dfs); | |
| 461 | } | ||
| 462 | ✗ | out[small_store].dofs.push_back(dfs); | |
| 463 | } | ||
| 464 | // Now we set the larger storage location to empty. | ||
| 465 | ✗ | out[large_store].dofs = {}; | |
| 466 | ✗ | out[large_store].coef = 0; | |
| 467 | } | ||
| 468 | } | ||
| 469 | } else { | ||
| 470 | // With the exception of corners, this will be the default case. | ||
| 471 | // We just add a pair of dofs to the bookkeeping. | ||
| 472 |
1/2✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
|
1360 | d.dofs.push_back(small_dof); |
| 473 |
1/2✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
|
1360 | d.dofs.push_back(large_dof); |
| 474 | 1360 | already_stored_in[small_dof] = d_count; | |
| 475 | 1360 | already_stored_in[large_dof] = d_count; | |
| 476 | 1360 | d_count++; | |
| 477 | 1360 | d.coef = 1; | |
| 478 |
1/2✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
|
1360 | out.push_back(d); |
| 479 | } | ||
| 480 | } | ||
| 481 | } | ||
| 482 | } | ||
| 483 | // Now we need to clean up the empty dofsets, since subsequent routines | ||
| 484 | // assume there to be at least one element. | ||
| 485 |
2/2✓ Branch 4 taken 1360 times.
✓ Branch 5 taken 66 times.
|
1426 | for (auto x = out.begin(); x != out.end(); ++x) { |
| 486 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 1360 times.
|
1360 | if ((*x).dofs.size() == 0) { |
| 487 | ✗ | out.erase(x); | |
| 488 | } | ||
| 489 | } | ||
| 490 | 132 | return out; | |
| 491 | 66 | } | |
| 492 | }; | ||
| 493 | |||
| 494 | /** | ||
| 495 | * \brief Specification of the struct for the vector Div conforming case. | ||
| 496 | * | ||
| 497 | * The Maxwell case: Here, the identification is 1-to-1, but we need to take | ||
| 498 | * care to glue the right vector component. The function builds upon the fact, | ||
| 499 | * that px == py == p. | ||
| 500 | */ | ||
| 501 | template <typename Derived> | ||
| 502 | struct glue_identificationmaker_<Derived, DifferentialForm::DivConforming> { | ||
| 503 | 67 | static std::vector<dofIdentification> makeIdentification( | |
| 504 | const std::vector<std::array<int, 4>> &edges_, | ||
| 505 | const SuperSpace<Derived> &superspace, const Projector<Derived> &proj) { | ||
| 506 | // since we assume px == py and uniform refinement in both directions, there | ||
| 507 | // will be a small_dim and a large_dim in every vector component. | ||
| 508 | 67 | const int small_dim = superspace.get_polynomial_degree() + | |
| 509 | 67 | proj.get_knot_repetition() * | |
| 510 | 67 | ((1 << superspace.get_refinement_level()) - 1); | |
| 511 | 67 | const int large_dim = superspace.get_polynomial_degree() + 1 + | |
| 512 | 67 | proj.get_knot_repetition() * | |
| 513 | 67 | ((1 << superspace.get_refinement_level()) - 1); | |
| 514 | 67 | const int dofs_per_patch_per_component = small_dim * large_dim; | |
| 515 | // sanity check | ||
| 516 |
1/2✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
|
67 | assert(dofs_per_patch_per_component == |
| 517 | (proj.get_dofs_after_projector() / | ||
| 518 | (superspace.get_number_of_patches() * 2)) && | ||
| 519 | "The assembly of the glue matrix is highly specific to the space. " | ||
| 520 | "Something went wrong; is the discrete space correct?"); | ||
| 521 | |||
| 522 | 67 | std::vector<GlueRoutines::dofIdentification> out; | |
| 523 |
1/2✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
|
67 | out.reserve(edges_.size() * small_dim); |
| 524 | |||
| 525 |
2/2✓ Branch 7 taken 1052 times.
✓ Branch 8 taken 67 times.
|
1119 | for (auto edge : edges_) { |
| 526 |
3/4✓ Branch 2 taken 772 times.
✓ Branch 3 taken 280 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 772 times.
|
1052 | assert(edge[0] <= edge[1] || edge[1] == -1); |
| 527 | // The DOFs that need to be glued on edges 0 and 2 are in the second | ||
| 528 | // vector component. The remainder of the dof-index-shift is to account | ||
| 529 | // for dofs on other patches. | ||
| 530 | 1052 | const int shift_e1 = edge[0] * dofs_per_patch_per_component + | |
| 531 | 1052 | (edgeToBeGluedInFirstComp(edge[2]) | |
| 532 |
2/2✓ Branch 0 taken 528 times.
✓ Branch 1 taken 524 times.
|
1576 | ? 0 |
| 533 | : (dofs_per_patch_per_component * | ||
| 534 | 524 | superspace.get_number_of_patches())); | |
| 535 | 1052 | const int shift_e2 = edge[1] * dofs_per_patch_per_component + | |
| 536 | 1052 | (edgeToBeGluedInFirstComp(edge[3]) | |
| 537 |
2/2✓ Branch 0 taken 910 times.
✓ Branch 1 taken 142 times.
|
1194 | ? 0 |
| 538 | : (dofs_per_patch_per_component * | ||
| 539 | 142 | superspace.get_number_of_patches())); | |
| 540 | |||
| 541 | 1052 | const int xdir_dim_1 = | |
| 542 |
2/2✓ Branch 2 taken 528 times.
✓ Branch 3 taken 524 times.
|
1052 | edgeToBeGluedInFirstComp(edge[2]) ? large_dim : small_dim; |
| 543 | 1052 | const int ydir_dim_1 = | |
| 544 |
2/2✓ Branch 2 taken 528 times.
✓ Branch 3 taken 524 times.
|
1052 | edgeToBeGluedInFirstComp(edge[2]) ? small_dim : large_dim; |
| 545 | 1052 | const int xdir_dim_2 = | |
| 546 |
2/2✓ Branch 2 taken 910 times.
✓ Branch 3 taken 142 times.
|
1052 | edgeToBeGluedInFirstComp(edge[3]) ? large_dim : small_dim; |
| 547 | 1052 | const int ydir_dim_2 = | |
| 548 |
2/2✓ Branch 2 taken 910 times.
✓ Branch 3 taken 142 times.
|
1052 | edgeToBeGluedInFirstComp(edge[3]) ? small_dim : large_dim; |
| 549 |
1/2✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
|
1052 | std::vector<int> dofs_e1 = |
| 550 | 1052 | getEdgeDofIndices(edge[2], xdir_dim_1, ydir_dim_1, shift_e1); | |
| 551 |
1/2✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
|
1052 | std::vector<int> dofs_e2 = |
| 552 | 1052 | getEdgeDofIndices(edge[3], xdir_dim_2, ydir_dim_2, shift_e2); | |
| 553 | |||
| 554 | 1052 | const int size_of_edge_dofs = dofs_e1.size(); | |
| 555 |
3/4✓ Branch 1 taken 772 times.
✓ Branch 2 taken 280 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 772 times.
|
1052 | assert(size_of_edge_dofs == dofs_e2.size() || edge[1] == -1); |
| 556 | |||
| 557 | 1052 | const bool needReversion = reverseParametrized(edge); | |
| 558 | 1052 | const int coef = glueCoefficientDivergenceConforming(edge); | |
| 559 | |||
| 560 | // Check if the edge is hanging. For screens one would need an "else". | ||
| 561 |
5/6✓ Branch 1 taken 280 times.
✓ Branch 2 taken 772 times.
✓ Branch 4 taken 280 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 280 times.
✓ Branch 7 taken 772 times.
|
1052 | if (edge[1] > -1 && edge[0] > -1) { |
| 562 |
2/2✓ Branch 1 taken 1404 times.
✓ Branch 2 taken 280 times.
|
1684 | for (int i = 0; i < size_of_edge_dofs; ++i) { |
| 563 | 1404 | GlueRoutines::dofIdentification d; | |
| 564 |
2/2✓ Branch 0 taken 702 times.
✓ Branch 1 taken 702 times.
|
1404 | const int j = needReversion ? size_of_edge_dofs - 1 - i : i; |
| 565 |
1/2✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
|
1404 | assert(dofs_e1[i] != dofs_e2[j] && |
| 566 | "If this happens something went horribly wrong."); | ||
| 567 |
1/2✓ Branch 4 taken 1404 times.
✗ Branch 5 not taken.
|
1404 | d.dofs.push_back(std::min(dofs_e1[i], dofs_e2[j])); |
| 568 |
1/2✓ Branch 4 taken 1404 times.
✗ Branch 5 not taken.
|
1404 | d.dofs.push_back(std::max(dofs_e1[i], dofs_e2[j])); |
| 569 | 1404 | d.coef = coef; | |
| 570 |
1/2✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
|
1404 | out.push_back(d); |
| 571 | } | ||
| 572 | } | ||
| 573 | } | ||
| 574 |
1/2✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
|
67 | out.shrink_to_fit(); |
| 575 | 67 | return out; | |
| 576 | ✗ | } | |
| 577 | }; | ||
| 578 | }; // namespace GlueRoutines | ||
| 579 | |||
| 580 | } // namespace Bembel | ||
| 581 | #endif // BEMBEL_SRC_ANSATZSPACE_GLUE_HPP_ | ||
| 582 |