GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/common/reorderingdofmapper.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 36 36 100.0%
Functions: 3 3 100.0%
Branches: 89 259 34.4%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7
8 /*!
9 * \file
10 * \ingroup Core
11 * \brief An SCSG element mapper that sorts the indices in order to optimize the matrix sparsity pattern
12 * \note The reordering needs the SCOTCH library
13 */
14 #ifndef DUMUX_COMMON_REORDERING_DOF_MAPPER_HH
15 #define DUMUX_COMMON_REORDERING_DOF_MAPPER_HH
16
17 #if HAVE_PTSCOTCH
18
19 #include <dune/common/timer.hh>
20 #include <dune/grid/common/mapper.hh>
21
22 #include <dumux/linear/scotchbackend.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup Core
28 * \brief An SCSG element mapper that sorts the indices in order to optimize the matrix sparsity pattern
29 * \note The reordering needs the SCOTCH library
30 */
31 template<class GridView>
32 7 class ReorderingDofMapper
33 : public Dune::Mapper<typename GridView::Grid, ReorderingDofMapper<GridView>, typename GridView::IndexSet::IndexType>
34 {
35 using Index = typename GridView::IndexSet::IndexType;
36 using ParentType = typename Dune::Mapper<typename GridView::Grid, ReorderingDofMapper<GridView>, Index>;
37 using Element = typename GridView::template Codim<0>::Entity;
38 public:
39
40 /*!
41 * \brief Construct mapper from grid and one of its index sets.
42 * \param gridView A Dune GridView object.
43 * \param layout a layout object (we just check whether it contains elements -> element mapper, else it's a vertex mapper)
44 */
45 template<class Layout>
46 13 ReorderingDofMapper (const GridView& gridView, Layout&& layout)
47
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 : gridView_(gridView)
48 13 , indexSet_(&gridView.indexSet())
49
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
✓ Branch 5 taken 13 times.
✗ Branch 6 not taken.
 26 , codimension_(layout(indexSet_->types(0)[0], GridView::dimension) ? 0 : GridView::dimension)
50 {
51
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 update_();
52 13 }
53
54 /*!
55 * \brief Map entity to array index.
56 *
57 * \tparam EntityType
58 * \param e Reference to codim \a EntityType entity.
59 * \return An index in the range 0 ... Max number of entities in set - 1.
60 */
61 template<class EntityType>
62 3595196 Index index (const EntityType& e) const
63 {
64 // map the index using the permutation obtained from the reordering algorithm
65
30/62
✗ Branch 2 not taken.
✓ Branch 3 taken 4032 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 8 taken 177251 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 17 taken 5827 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 3751 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 1720 times.
✓ Branch 23 taken 4107 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 2016 times.
✓ Branch 27 taken 184687 times.
✓ Branch 31 taken 20344 times.
✓ Branch 32 taken 1980 times.
✓ Branch 34 taken 3470 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 7348 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 14010 times.
✓ Branch 40 taken 1 times.
✓ Branch 41 taken 1728 times.
✓ Branch 42 taken 6055 times.
✓ Branch 44 taken 1735 times.
✗ Branch 45 not taken.
✓ Branch 48 taken 1735 times.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✓ Branch 51 taken 3674 times.
✓ Branch 53 taken 1735 times.
✗ Branch 54 not taken.
✗ Branch 59 not taken.
✓ Branch 60 taken 3470 times.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✓ Branch 66 taken 2873 times.
✗ Branch 67 not taken.
✓ Branch 69 taken 1735 times.
✗ Branch 70 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✓ Branch 76 taken 1735 times.
✗ Branch 77 not taken.
✗ Branch 79 not taken.
✓ Branch 80 taken 2873 times.
✓ Branch 83 taken 1735 times.
✗ Branch 84 not taken.
✗ Branch 4 not taken.
✗ Branch 10 not taken.
✓ Branch 13 taken 2016 times.
✗ Branch 14 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 4092 times.
✓ Branch 30 taken 12276 times.
✓ Branch 33 taken 36 times.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
 2173812 return static_cast<Index>(permutation_[indexSet_->index(e)]);
66 }
67
68 /** @brief Map subentity of codim 0 entity to array index.
69
70 \param e Reference to codim 0 entity.
71 \param i Number of subentity of e
72 \param codim Codimension of the subentity
73 \return An index in the range 0 ... Max number of entities in set - 1.
74 */
75 84662 Index subIndex (const Element& e, int i, unsigned int codim) const
76 {
77
1/8
✓ Branch 6 taken 4032 times.
✗ Branch 7 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 84662 return indexSet_->subIndex(e, i, codim);
78 }
79
80 /** @brief Return total number of entities in the entity set managed by the mapper.
81
82 This number can be used to allocate a vector of data elements associated with the
83 entities of the set. In the parallel case this number is per process (i.e. it
84 may be different in different processes).
85
86 \return Size of the entity set.
87 */
88 760 std::size_t size () const
89 {
90
20/128
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 24 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✓ Branch 40 taken 24 times.
✓ Branch 41 taken 101 times.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 101 times.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✓ Branch 62 taken 12 times.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✓ Branch 65 taken 12 times.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✓ Branch 72 taken 12 times.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✓ Branch 75 taken 24 times.
✗ Branch 76 not taken.
✓ Branch 77 taken 12 times.
✗ Branch 78 not taken.
✓ Branch 80 taken 6 times.
✗ Branch 81 not taken.
✓ Branch 82 taken 508 times.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✓ Branch 86 taken 508 times.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✓ Branch 94 taken 6 times.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
✓ Branch 98 taken 6 times.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✓ Branch 102 taken 6 times.
✗ Branch 103 not taken.
✗ Branch 105 not taken.
✓ Branch 106 taken 6 times.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✓ Branch 109 taken 6 times.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 114 not taken.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
✗ Branch 120 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✓ Branch 127 taken 3 times.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✓ Branch 130 taken 3 times.
✗ Branch 131 not taken.
✗ Branch 54 not taken.
✗ Branch 57 not taken.
✗ Branch 68 not taken.
✗ Branch 79 not taken.
✓ Branch 90 taken 6 times.
✗ Branch 93 not taken.
✗ Branch 104 not taken.
✗ Branch 115 not taken.
✗ Branch 121 not taken.
 1416 return indexSet_->size(codimension_);
91 }
92
93 /** @brief Returns true if the entity is contained in the index set
94
95 \param e Reference to entity
96 \param result integer reference where corresponding index is stored if true
97 \return true if entity is in entity set of the mapper
98 */
99 template<class EntityType>
100 bool contains (const EntityType& e, Index& result) const
101 {
102 result = index(e);
103 return true;
104 }
105
106 /** @brief Returns true if the entity is contained in the index set
107
108 \param e Reference to codim 0 entity
109 \param i subentity number
110 \param cc subentity codim
111 \param result integer reference where corresponding index is stored if true
112 \return true if entity is in entity set of the mapper
113 */
114 bool contains (const Element& e, int i, int cc, Index& result) const
115 {
116 result = indexSet_->subIndex(e, i, cc);
117 return true;
118 }
119
120 /*!
121 * \brief Recalculates map after mesh adaptation
122 */
123 13 void update (const GridView& gridView)
124 {
125
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 gridView_ = gridView;
126
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 indexSet_ = &gridView_.indexSet();
127 13 update_();
128 13 }
129
130 void update (GridView&& gridView)
131 {
132 gridView_ = std::move(gridView);
133 indexSet_ = &gridView_.indexSet();
134 update_();
135 }
136
137 private:
138 26 void update_()
139 {
140 // Compute scotch reordering
141
2/3
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 26 Dune::Timer watch;
142 // Create the graph as an adjacency list
143
2/3
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 52 std::vector<std::vector<Index>> graph(size());
144
145 // dofs on element centers (cell-centered methods)
146
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 12 times.
 26 if (codimension_ == 0)
147 {
148
2/2
✓ Branch 1 taken 7502 times.
✓ Branch 2 taken 14 times.
 7516 for (const auto& element : elements(gridView_))
149 {
150
1/2
✓ Branch 1 taken 7502 times.
✗ Branch 2 not taken.
 7502 auto eIdx = indexSet_->index(element);
151
5/8
✓ Branch 1 taken 7502 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 15844 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 23346 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 15532 times.
✓ Branch 11 taken 312 times.
 77228 for (const auto& intersection : intersections(gridView_, element))
152 {
153
2/2
✓ Branch 0 taken 15532 times.
✓ Branch 1 taken 312 times.
 15844 if (intersection.neighbor())
154
1/2
✓ Branch 1 taken 15532 times.
✗ Branch 2 not taken.
 15532 graph[eIdx].push_back(indexSet_->index(intersection.outside()));
155 }
156 }
157 }
158 // dof on vertices (box method)
159 else
160 {
161
2/2
✓ Branch 2 taken 4032 times.
✓ Branch 3 taken 12 times.
 8076 for (const auto& element : elements(gridView_))
162 {
163 4032 auto eIdx = indexSet_->index(element);
164
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 12096 times.
✓ Branch 2 taken 8064 times.
✓ Branch 3 taken 4032 times.
 24192 for (int vIdxLocal = 0; vIdxLocal < element.subEntities(codimension_); ++vIdxLocal)
165 {
166 8064 auto vIdxGlobal = indexSet_->subIndex(element, vIdxLocal, codimension_);
167
1/2
✓ Branch 1 taken 8064 times.
✗ Branch 2 not taken.
 8064 graph[vIdxGlobal].push_back(eIdx);
168 }
169 }
170 }
171
172
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
 26 permutation_ = ScotchBackend<Index>::computeGPSReordering(graph);
173
4/7
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
 52 std::cout << "Scotch backend reordered index set of size " << size()
174
3/6
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
 26 << " in " << watch.elapsed() << " seconds." << std::endl;
175 26 }
176 // GridView is needed to keep the IndexSet valid
177 GridView gridView_;
178 const typename GridView::IndexSet* indexSet_;
179 const int codimension_;
180 // the map resulting from the reordering
181 std::vector<int> permutation_;
182 };
183
184 } // end namespace Dumux
185
186 #else
187
188 #warning "PTSCOTCH was not found on your system. Dumux::ReorderingDofMapper needs it to work -> fallback to MCMGMapper without reordering!"
189 #include <dune/grid/common/mcmgmapper.hh>
190 namespace Dumux {
191 template<class GridView>
192 using ReorderingDofMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
193 } // end namespace Dumux
194
195 #endif // HAVE_PTSCOTCH
196 #endif
197