Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #ifndef LIBMESH_SPARSITY_PATTERN_H
20 : #define LIBMESH_SPARSITY_PATTERN_H
21 :
22 : // Local Includes
23 : #include "libmesh/elem_range.h"
24 : #include "libmesh/threads_allocators.h"
25 : #include "libmesh/parallel_object.h"
26 :
27 : // C++ includes
28 : #include <algorithm> // is_sorted
29 : #include <unordered_set>
30 : #include <vector>
31 :
32 : namespace libMesh
33 : {
34 :
35 : // Forward declarations
36 : class DofMap;
37 : class CouplingMatrix;
38 : class StaticCondensationDofMap;
39 :
40 : /**
41 : * This defines the sparsity pattern, or graph, of a sparse matrix.
42 : * The format is quite simple -- the global indices of the nonzero entries
43 : * in each row are packed into a vector. The global indices (i,j) of the
44 : * nth nonzero entry of row i are given by j = sparsity_pattern[i][n];
45 : *
46 : * \author Roy Stogner
47 : * \date 2010
48 : * \brief Defines the sparsity pattern of a sparse matrix.
49 : */
50 : namespace SparsityPattern // use a namespace so member classes can be forward-declared.
51 : {
52 : typedef std::vector<dof_id_type, Threads::scalable_allocator<dof_id_type>> Row;
53 39481 : class Graph : public std::vector<Row> {};
54 :
55 39481 : class NonlocalGraph : public std::map<dof_id_type, Row> {};
56 :
57 : /**
58 : * Splices the two sorted ranges [begin,middle) and [middle,end)
59 : * into one sorted range [begin,end). This method is much like
60 : * std::inplace_merge except it assumes the intersection
61 : * of the two sorted ranges is empty and that any element in
62 : * each range occurs only once in that range. Additionally,
63 : * this sort occurs in-place, while std::inplace_merge may
64 : * use a temporary buffer.
65 : */
66 : template<typename BidirectionalIterator>
67 : static void sort_row (const BidirectionalIterator begin,
68 : BidirectionalIterator middle,
69 : const BidirectionalIterator end);
70 :
71 :
72 :
73 : /**
74 : * Abstract base class to be used to add user-defined implicit
75 : * degree of freedom couplings.
76 : */
77 : class AugmentSparsityPattern
78 : {
79 : public:
80 : virtual ~AugmentSparsityPattern () = default;
81 :
82 : /**
83 : * User-defined function to augment the sparsity pattern.
84 : */
85 : virtual void augment_sparsity_pattern (SparsityPattern::Graph & sparsity,
86 : std::vector<dof_id_type> & n_nz,
87 : std::vector<dof_id_type> & n_oz) = 0;
88 : };
89 :
90 :
91 :
92 : /**
93 : * This helper class can be called on multiple threads to compute
94 : * the sparsity pattern (or graph) of the sparse matrix resulting
95 : * from the discretization. This pattern may be used directly by
96 : * a particular sparse matrix format (e.g. \p LaspackMatrix)
97 : * or indirectly (e.g. \p PetscMatrixBase). In the latter case the
98 : * number of nonzeros per row of the matrix is needed for efficient
99 : * preallocation. In this case it suffices to provide estimate
100 : * (but bounding) values, and in this case the threaded method can
101 : * take some short-cuts for efficiency.
102 : */
103 : class Build : public ParallelObject
104 : {
105 : public:
106 : Build (const DofMap & dof_map_in,
107 : const CouplingMatrix * dof_coupling_in,
108 : const std::vector<GhostingFunctor *> & coupling_functors_in,
109 : const bool implicit_neighbor_dofs_in,
110 : const bool need_full_sparsity_pattern_in,
111 : const bool calculate_constrained_in = false,
112 : const StaticCondensationDofMap * sc = nullptr);
113 :
114 : /**
115 : * Special functions.
116 : * - The Build object holds references to DofMap and coupling
117 : * functors, therefore it can't be assigned.
118 : */
119 : Build (const Build &) = default;
120 : Build & operator= (const Build &) = delete;
121 : Build (Build &&) = default;
122 : Build & operator= (Build &&) = delete;
123 83358 : ~Build () = default;
124 :
125 : /**
126 : * Splitting constructor, for use in multithreaded loops.
127 : */
128 : Build (Build & other, Threads::split);
129 :
130 : /**
131 : * Add entries from a range of elements to this object's sparsity
132 : * pattern.
133 : */
134 : void operator()(const ConstElemRange & range);
135 :
136 : /**
137 : * Combine the sparsity pattern in \p other with this object's
138 : * sparsity pattern. Useful in multithreaded loops.
139 : */
140 : void join (const Build & other);
141 :
142 : /**
143 : * Send sparsity pattern data relevant to other processors to those
144 : * processors, and receive and incorporate data relevant to us.
145 : */
146 : void parallel_sync ();
147 :
148 : /**
149 : * Rows of sparse matrix indices, indexed by the offset from the
150 : * first DoF on this processor.
151 : */
152 1240 : const SparsityPattern::Graph & get_sparsity_pattern() const
153 1240 : { return sparsity_pattern; }
154 :
155 : /**
156 : * Rows of sparse matrix indices, mapped from global DoF number,
157 : * which belong on other processors. Stored here only temporarily
158 : * until a parallel_sync() sends them where they belong.
159 : */
160 : const SparsityPattern::NonlocalGraph & get_nonlocal_pattern() const
161 : { return nonlocal_pattern; }
162 :
163 : /**
164 : * The total number of nonzeros in the global matrix.
165 : */
166 : std::size_t n_nonzeros() const;
167 :
168 : /**
169 : * The number of on-processor nonzeros in my portion of the
170 : * global matrix.
171 : */
172 2270 : const std::vector<dof_id_type> & get_n_nz() const
173 43634 : { return n_nz; }
174 :
175 : /**
176 : * The number of off-processor nonzeros in my portion of the
177 : * global matrix.
178 : */
179 2218 : const std::vector<dof_id_type> & get_n_oz() const
180 41370 : { return n_oz; }
181 :
182 : /**
183 : * Let a user-provided AugmentSparsityPattern subclass modify our
184 : * sparsity structure.
185 : */
186 : void apply_extra_sparsity_object(SparsityPattern::AugmentSparsityPattern & asp);
187 :
188 : /**
189 : * Let a user-provided function modify our sparsity structure.
190 : */
191 0 : void apply_extra_sparsity_function(void (*func)(SparsityPattern::Graph & sparsity,
192 : std::vector<dof_id_type> & n_nz,
193 : std::vector<dof_id_type> & n_oz,
194 : void * context),
195 : void * context)
196 0 : { func(sparsity_pattern, n_nz, n_oz, context); }
197 :
198 : /**
199 : * Clear the "full" details of our sparsity structure, leaving only
200 : * the counts of non-zero entries.
201 : */
202 37861 : void clear_full_sparsity()
203 : {
204 36651 : sparsity_pattern.clear();
205 1210 : nonlocal_pattern.clear();
206 37861 : }
207 :
208 : private:
209 : const DofMap & dof_map;
210 : const CouplingMatrix * dof_coupling;
211 : const std::vector<GhostingFunctor *> & coupling_functors;
212 : const bool implicit_neighbor_dofs;
213 : const bool need_full_sparsity_pattern;
214 : const bool calculate_constrained;
215 : const StaticCondensationDofMap * const sc;
216 :
217 : /**
218 : * If there are "spider" nodes in the mesh (i.e. a single node which
219 : * is connected to many 1D elements) and Constraints, we can end up
220 : * sorting the same set of DOFs multiple times in handle_vi_vj(),
221 : * only to find that it has no net effect on the final sparsity. In
222 : * such cases it is much faster to keep track of (element_dofs_i,
223 : * element_dofs_j) pairs which have already been handled and not
224 : * repeat the computation. We use this data structure to keep track
225 : * of hashes of sets of dofs we have already seen, thus avoiding
226 : * unnecessary calculations.
227 : */
228 : std::unordered_set<dof_id_type> hashed_dof_sets;
229 :
230 : /// A dummy work vector to avoid repeated memory allocations
231 : std::vector<dof_id_type> dummy_vec;
232 :
233 : void handle_vi_vj(const std::vector<dof_id_type> & element_dofs_i,
234 : const std::vector<dof_id_type> & element_dofs_j);
235 :
236 : void sorted_connected_dofs(const Elem * elem,
237 : std::vector<dof_id_type> & dofs_vi,
238 : unsigned int vi);
239 :
240 : #ifndef LIBMESH_ENABLE_DEPRECATED
241 : private:
242 : #endif
243 :
244 : SparsityPattern::Graph sparsity_pattern;
245 :
246 : SparsityPattern::NonlocalGraph nonlocal_pattern;
247 :
248 : std::vector<dof_id_type> n_nz;
249 :
250 : std::vector<dof_id_type> n_oz;
251 : };
252 :
253 : }
254 :
255 :
256 :
257 : // ------------------------------------------------------------
258 : // SparsityPattern inline member functions
259 : template<typename BidirectionalIterator>
260 : inline
261 39346947 : void SparsityPattern::sort_row (const BidirectionalIterator begin,
262 : BidirectionalIterator middle,
263 : const BidirectionalIterator end)
264 : {
265 : // Assure we have the conditions for an inplace_merge
266 : #ifdef DEBUG
267 3973576 : libmesh_assert(std::is_sorted(begin, middle));
268 3973576 : libmesh_assert(std::is_sorted(middle, end));
269 : #endif
270 3973576 : libmesh_assert(std::unique(begin, middle) == middle);
271 3973576 : libmesh_assert(std::unique(middle, end) == end);
272 :
273 43320523 : std::inplace_merge(begin, middle, end);
274 :
275 : // Assure the algorithm worked if we are in DEBUG mode
276 : #ifdef DEBUG
277 3973576 : libmesh_assert (std::is_sorted(begin,end));
278 : #endif
279 :
280 : // Make sure the two ranges did not contain any common elements
281 3973576 : libmesh_assert (std::unique (begin, end) == end);
282 39346947 : }
283 :
284 : } // namespace libMesh
285 :
286 : #endif // LIBMESH_SPARSITY_PATTERN_H
|
