Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 :
20 : // Local includes
21 : #include "libmesh/sparsity_pattern.h"
22 :
23 : // libMesh includes
24 : #include "libmesh/coupling_matrix.h"
25 : #include "libmesh/dof_map.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/ghosting_functor.h"
28 : #include "libmesh/hashword.h"
29 : #include "libmesh/parallel_algebra.h"
30 : #include "libmesh/parallel.h"
31 : #include "libmesh/parallel_sync.h"
32 : #include "libmesh/utility.h"
33 : #include "libmesh/reduced_dof_map.h"
34 :
35 : // TIMPI includes
36 : #include "timpi/communicator.h"
37 :
38 :
39 : namespace libMesh
40 : {
41 : namespace SparsityPattern
42 : {
43 :
44 : //-------------------------------------------------------
45 : // we need to implement these constructors here so that
46 : // a full DofMap definition is available.
47 4680 : Build::Build (const DofMap & dof_map_in,
48 : const CouplingMatrix * dof_coupling_in,
49 : const std::vector<GhostingFunctor *> & coupling_functors_in,
50 : const bool implicit_neighbor_dofs_in,
51 : const bool need_full_sparsity_pattern_in,
52 : const bool calculate_constrained_in,
53 4680 : const ReducedDofMap * const reduced_in) :
54 : ParallelObject(dof_map_in),
55 1852 : dof_map(dof_map_in),
56 1852 : dof_coupling(dof_coupling_in),
57 1852 : coupling_functors(coupling_functors_in),
58 1852 : implicit_neighbor_dofs(implicit_neighbor_dofs_in),
59 1852 : need_full_sparsity_pattern(need_full_sparsity_pattern_in),
60 1852 : calculate_constrained(calculate_constrained_in),
61 1852 : reduced(reduced_in),
62 : sparsity_pattern(),
63 : nonlocal_pattern(),
64 : n_nz(),
65 4680 : n_oz()
66 4680 : {}
67 :
68 :
69 :
70 84 : Build::Build (Build & other, Threads::split) :
71 : ParallelObject(other),
72 84 : dof_map(other.dof_map),
73 84 : dof_coupling(other.dof_coupling),
74 84 : coupling_functors(other.coupling_functors),
75 84 : implicit_neighbor_dofs(other.implicit_neighbor_dofs),
76 84 : need_full_sparsity_pattern(other.need_full_sparsity_pattern),
77 84 : calculate_constrained(other.calculate_constrained),
78 84 : reduced(other.reduced),
79 28 : hashed_dof_sets(other.hashed_dof_sets),
80 : sparsity_pattern(),
81 : nonlocal_pattern(),
82 : n_nz(),
83 84 : n_oz()
84 84 : {}
85 :
86 :
87 :
88 : #if defined(__GNUC__) && (__GNUC__ < 4) && !defined(__INTEL_COMPILER)
89 :
90 : void _dummy_function(void) {}
91 :
92 : #endif
93 :
94 :
95 :
96 2732265 : void Build::sorted_connected_dofs(const Elem * elem,
97 : std::vector<dof_id_type> & dofs_vi,
98 : unsigned int vi)
99 : {
100 2732265 : if (this->reduced)
101 : {
102 : // We build a sparsity pattern that will match the size of the reduced system. This is so that
103 : // we have the data necessary to init the reduced system matrix
104 47061 : dofs_vi.clear();
105 :
106 : // Unretained dofs still couple their constraining dofs, which constrained element matrices
107 : // and vectors act on
108 47061 : dummy_unretained.clear();
109 :
110 : auto total_and_retained_from_scalar_dofs_functor =
111 16 : [&dofs_vi, this](const Elem & /*elem*/,
112 : std::vector<dof_id_type> & dof_indices,
113 80 : const std::vector<dof_id_type> & scalar_dof_indices)
114 : {
115 64 : dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
116 96 : for (const auto scalar_dof : scalar_dof_indices)
117 48 : (this->reduced->retains_scalar_dof(scalar_dof) ? dofs_vi : dummy_unretained)
118 48 : .push_back(scalar_dof);
119 94170 : };
120 :
121 : auto total_and_retained_from_field_dofs_functor =
122 612767 : [&dofs_vi, this](const Elem & functor_elem,
123 : const unsigned int node_num,
124 : const unsigned int var_num,
125 : std::vector<dof_id_type> & dof_indices,
126 2691744 : const dof_id_type field_dof)
127 : {
128 1838301 : dof_indices.push_back(field_dof);
129 1838301 : (this->reduced->retains_field_dof(functor_elem, node_num, var_num, field_dof)
130 : ? dofs_vi
131 : : dummy_unretained)
132 1838301 : .push_back(field_dof);
133 1979484 : };
134 :
135 329427 : dof_map.dof_indices(elem,
136 141183 : dummy_vec,
137 : vi,
138 : total_and_retained_from_scalar_dofs_functor,
139 : total_and_retained_from_field_dofs_functor,
140 94122 : elem->p_level());
141 :
142 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
143 : // Constraint rows reference unconstrained dofs only, so one level suffices
144 141183 : const auto & constraints = dof_map.get_dof_constraints();
145 502197 : for (const auto unretained : dummy_unretained)
146 361014 : if (const auto it = constraints.find(unretained); it != constraints.end())
147 2025 : for (const auto & [constraining_dof, weight] : it->second)
148 : {
149 360 : libmesh_ignore(weight);
150 1080 : dofs_vi.push_back(constraining_dof);
151 : }
152 : #endif
153 : }
154 : else
155 2591082 : dof_map.dof_indices (elem, dofs_vi, vi);
156 :
157 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
158 2732265 : dof_map.find_connected_dofs (dofs_vi);
159 : #endif
160 : // We can be more efficient if we sort the element DOFs into
161 : // increasing order
162 2732265 : std::sort(dofs_vi.begin(), dofs_vi.end());
163 :
164 : // Handle cases where duplicate nodes are intentionally assigned to
165 : // a single element.
166 2732265 : dofs_vi.erase(std::unique(dofs_vi.begin(), dofs_vi.end()), dofs_vi.end());
167 2732265 : }
168 :
169 :
170 :
171 4152577 : void Build::handle_vi_vj(const std::vector<dof_id_type> & element_dofs_i,
172 : const std::vector<dof_id_type> & element_dofs_j)
173 : {
174 : const unsigned int n_dofs_on_element_i =
175 2237512 : cast_int<unsigned int>(element_dofs_i.size());
176 :
177 4152577 : const processor_id_type proc_id = dof_map.processor_id();
178 1118756 : const dof_id_type first_dof_on_proc = dof_map.first_dof(proc_id);
179 1118756 : const dof_id_type end_dof_on_proc = dof_map.end_dof(proc_id);
180 :
181 : std::vector<dof_id_type>
182 2237512 : dofs_to_add;
183 :
184 : const unsigned int n_dofs_on_element_j =
185 2237512 : cast_int<unsigned int>(element_dofs_j.size());
186 :
187 : // It only makes sense to compute hashes and see if we can skip
188 : // doing work when there are a "large" amount of DOFs for a given
189 : // element. The cutoff for "large" is somewhat arbitrarily chosen
190 : // based on a test case with a spider node that resulted in O(10^3)
191 : // entries in element_dofs_i for O(10^3) elements. Making this
192 : // number larger will disable the hashing optimization in more
193 : // cases.
194 1118756 : bool dofs_seen = false;
195 4152577 : if (n_dofs_on_element_j > 0 && n_dofs_on_element_i > 256)
196 : {
197 0 : auto hash_i = Utility::hashword(element_dofs_i);
198 0 : auto hash_j = Utility::hashword(element_dofs_j);
199 0 : auto final_hash = Utility::hashword2(hash_i, hash_j);
200 0 : auto result = hashed_dof_sets.insert(final_hash);
201 : // if insert failed, we have already seen these dofs
202 0 : dofs_seen = !result.second;
203 : }
204 :
205 : // there might be 0 dofs for the other variable on the same element
206 : // (when subdomain variables do not overlap) and that's when we do
207 : // not do anything
208 4152577 : if (n_dofs_on_element_j > 0 && !dofs_seen)
209 : {
210 27249623 : for (unsigned int i=0; i<n_dofs_on_element_i; i++)
211 : {
212 23114332 : const dof_id_type ig = element_dofs_i[i];
213 :
214 : SparsityPattern::Row * row;
215 :
216 : // We save non-local row components for now so we can
217 : // communicate them to other processors later.
218 :
219 23114332 : if ((ig >= first_dof_on_proc) &&
220 6101633 : (ig < end_dof_on_proc))
221 : {
222 : // This is what I mean
223 : // libmesh_assert_greater_equal ((ig - first_dof_on_proc), 0);
224 : // but do the test like this because ig and
225 : // first_dof_on_proc are unsigned ints
226 5972553 : libmesh_assert_greater_equal (ig, first_dof_on_proc);
227 5972553 : libmesh_assert_less (ig, (sparsity_pattern.size() +
228 : first_dof_on_proc));
229 :
230 28309673 : row = &sparsity_pattern[ig - first_dof_on_proc];
231 : }
232 : else
233 : {
234 777212 : row = &nonlocal_pattern[ig];
235 : }
236 :
237 : // If the row is empty we will add *all*
238 : // the element j DOFs, so just do that.
239 23114332 : if (row->empty())
240 : {
241 2343961 : row->insert(row->end(),
242 : element_dofs_j.begin(),
243 3445604 : element_dofs_j.end());
244 : }
245 : else
246 : {
247 : // Build a list of the DOF indices not found in the
248 : // sparsity pattern
249 5378178 : dofs_to_add.clear();
250 :
251 : // Cache iterators. Low will move forward, subsequent
252 : // searches will be on smaller ranges
253 : SparsityPattern::Row::iterator
254 : low = std::lower_bound
255 14530792 : (row->begin(), row->end(), element_dofs_j.front()),
256 : high = std::upper_bound
257 14530792 : (low, row->end(), element_dofs_j.back());
258 :
259 192350096 : for (unsigned int j=0; j<n_dofs_on_element_j; j++)
260 : {
261 218684470 : const dof_id_type jg = element_dofs_j[j];
262 :
263 : // See if jg is in the sorted range
264 : std::pair<SparsityPattern::Row::iterator,
265 : SparsityPattern::Row::iterator>
266 172441126 : pos = std::equal_range (low, high, jg);
267 :
268 : // Must add jg if it wasn't found
269 172441126 : if (pos.first == pos.second)
270 90513124 : dofs_to_add.push_back(jg);
271 :
272 : // pos.first is now a valid lower bound for any
273 : // remaining element j DOFs. (That's why we sorted them.)
274 : // Use it for the next search
275 46243344 : low = pos.first;
276 : }
277 :
278 : // Add to the sparsity pattern
279 19908970 : if (!dofs_to_add.empty())
280 : {
281 8640450 : const std::size_t old_size = row->size();
282 :
283 11843692 : row->insert (row->end(),
284 : dofs_to_add.begin(),
285 12960675 : dofs_to_add.end());
286 :
287 : SparsityPattern::sort_row
288 11843692 : (row->begin(), row->begin()+old_size,
289 : row->end());
290 : }
291 : }
292 : } // End dofs-of-var-i loop
293 : } // End if-dofs-of-var-j
294 4152577 : }
295 :
296 :
297 :
298 4764 : void Build::operator()(const ConstElemRange & range)
299 : {
300 : // Compute the sparsity structure of the global matrix. This can be
301 : // fed into a PetscMatrixBase to allocate exactly the number of nonzeros
302 : // necessary to store the matrix. This algorithm should be linear
303 : // in the (# of elements)*(# nodes per element)
304 6206 : sparsity_pattern.resize(dof_map.n_local_dofs());
305 :
306 : // Handle dof coupling specified by library and user coupling functors
307 : {
308 4764 : const unsigned int n_var = dof_map.n_variables();
309 :
310 7648 : std::vector<std::vector<dof_id_type> > element_dofs_i(n_var);
311 :
312 2884 : std::vector<const Elem *> coupled_neighbors;
313 1191323 : for (const auto & elem : range)
314 : {
315 : // Make some fake element iterators defining a range
316 : // pointing to only this element.
317 1186559 : Elem * const * elempp = const_cast<Elem * const *>(&elem);
318 1186559 : Elem * const * elemend = elempp+1;
319 :
320 : const MeshBase::const_element_iterator fake_elem_it =
321 : MeshBase::const_element_iterator(elempp,
322 : elemend,
323 1507654 : Predicates::NotNull<Elem * const *>());
324 :
325 : const MeshBase::const_element_iterator fake_elem_end =
326 : MeshBase::const_element_iterator(elemend,
327 : elemend,
328 2373118 : Predicates::NotNull<Elem * const *>());
329 :
330 642190 : GhostingFunctor::map_type elements_to_couple;
331 642190 : DofMap::CouplingMatricesSet temporary_coupling_matrices;
332 :
333 1186559 : dof_map.merge_ghost_functor_outputs(elements_to_couple,
334 : temporary_coupling_matrices,
335 1186559 : dof_map.coupling_functors_begin(),
336 1507654 : dof_map.coupling_functors_end(),
337 : fake_elem_it,
338 : fake_elem_end,
339 : DofObject::invalid_processor_id);
340 2776340 : for (unsigned int vi=0; vi<n_var; vi++)
341 2016927 : this->sorted_connected_dofs(elem, element_dofs_i[vi], vi);
342 :
343 2776340 : for (unsigned int vi=0; vi<n_var; vi++)
344 4308114 : for (const auto & [partner, ghost_coupling] : elements_to_couple)
345 : {
346 : // Loop over coupling matrix row variables if we have a
347 : // coupling matrix, or all variables if not.
348 2718333 : if (ghost_coupling)
349 : {
350 16 : libmesh_assert_equal_to (ghost_coupling->size(), n_var);
351 48 : ConstCouplingRow ccr(vi, *ghost_coupling);
352 :
353 108 : for (const auto & idx : ccr)
354 : {
355 60 : if (partner == elem)
356 0 : this->handle_vi_vj(element_dofs_i[vi], element_dofs_i[idx]);
357 : else
358 : {
359 40 : std::vector<dof_id_type> partner_dofs;
360 60 : this->sorted_connected_dofs(partner, partner_dofs, idx);
361 80 : this->handle_vi_vj(element_dofs_i[vi], partner_dofs);
362 : }
363 : }
364 : }
365 : else
366 : {
367 6870802 : for (unsigned int vj = 0; vj != n_var; ++vj)
368 : {
369 4152517 : if (partner == elem)
370 4608225 : this->handle_vi_vj(element_dofs_i[vi], element_dofs_i[vj]);
371 : else
372 : {
373 639340 : std::vector<dof_id_type> partner_dofs;
374 1142424 : this->sorted_connected_dofs(partner, partner_dofs, vj);
375 1462094 : this->handle_vi_vj(element_dofs_i[vi], partner_dofs);
376 : }
377 : }
378 : }
379 : } // End ghosted element loop
380 : } // End range element loop
381 1880 : } // End ghosting functor section
382 4764 : }
383 :
384 :
385 :
386 84 : void Build::join (const SparsityPattern::Build & other)
387 : {
388 28 : libmesh_assert_equal_to (sparsity_pattern.size(), other.sparsity_pattern.size());
389 :
390 323247 : for (dof_id_type r=0; r<dof_map.n_local_dofs(); r++)
391 : {
392 : // increment the number of on and off-processor nonzeros in this row
393 : // (note this will be an upper bound unless we need the full sparsity pattern)
394 211040 : SparsityPattern::Row & my_row = sparsity_pattern[r];
395 211040 : const SparsityPattern::Row & their_row = other.sparsity_pattern[r];
396 :
397 : // simple copy if I have no dofs
398 323135 : if (my_row.empty())
399 143876 : my_row = their_row;
400 :
401 : // otherwise add their DOFs to mine, resort, and re-unique the row
402 179259 : else if (!their_row.empty()) // do nothing for the trivial case where
403 : { // their row is empty
404 20390 : my_row.insert (my_row.end(),
405 : their_row.begin(),
406 30498 : their_row.end());
407 :
408 : // We cannot use SparsityPattern::sort_row() here because it expects
409 : // the [begin,middle) [middle,end) to be non-overlapping. This is not
410 : // necessarily the case here, so use std::sort()
411 30556 : std::sort (my_row.begin(), my_row.end());
412 :
413 30556 : my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
414 : }
415 : }
416 :
417 : // Move nonlocal row information to ourselves; the other thread
418 : // won't need it in the map after that.
419 3232 : for (const auto & p : other.nonlocal_pattern)
420 : {
421 : #ifndef NDEBUG
422 1043 : const dof_id_type dof_id = p.first;
423 :
424 1043 : processor_id_type dbg_proc_id = 0;
425 1056 : while (dof_id >= dof_map.end_dof(dbg_proc_id))
426 13 : dbg_proc_id++;
427 1043 : libmesh_assert (dbg_proc_id != this->processor_id());
428 : #endif
429 :
430 3148 : const SparsityPattern::Row & their_row = p.second;
431 :
432 : // We should have no empty values in a map
433 1043 : libmesh_assert (!their_row.empty());
434 :
435 3148 : if (auto my_it = nonlocal_pattern.find(p.first);
436 1043 : my_it == nonlocal_pattern.end())
437 : {
438 : // nonlocal_pattern[it->first].swap(their_row);
439 2686 : nonlocal_pattern[p.first] = their_row;
440 : }
441 : else
442 : {
443 462 : SparsityPattern::Row & my_row = my_it->second;
444 :
445 307 : my_row.insert (my_row.end(),
446 : their_row.begin(),
447 620 : their_row.end());
448 :
449 : // We cannot use SparsityPattern::sort_row() here because it expects
450 : // the [begin,middle) [middle,end) to be non-overlapping. This is not
451 : // necessarily the case here, so use std::sort()
452 462 : std::sort (my_row.begin(), my_row.end());
453 :
454 462 : my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
455 : }
456 : }
457 :
458 : // Combine the other thread's hashed_dof_sets with ours.
459 56 : hashed_dof_sets.insert(other.hashed_dof_sets.begin(),
460 : other.hashed_dof_sets.end());
461 84 : }
462 :
463 :
464 :
465 4680 : void Build::parallel_sync ()
466 : {
467 1414 : parallel_object_only();
468 1414 : libmesh_assert(this->comm().verify(need_full_sparsity_pattern));
469 :
470 4680 : const auto n_dofs_on_proc = dof_map.n_local_dofs();
471 1414 : const auto local_first_dof = dof_map.first_dof();
472 :
473 : // The data to send
474 2828 : std::map<processor_id_type, std::vector<dof_id_type>> ids_to_send;
475 2828 : std::map<processor_id_type, std::vector<Row>> rows_to_send;
476 :
477 : // Loop over the nonlocal rows and transform them into the new datastructure
478 1414 : NonlocalGraph::iterator it = nonlocal_pattern.begin();
479 93225 : while (it != nonlocal_pattern.end())
480 : {
481 88545 : const auto dof_id = it->first;
482 88545 : auto & row = it->second;
483 :
484 88545 : processor_id_type proc_id = 0;
485 139698 : while (dof_id >= dof_map.end_dof(proc_id))
486 15234 : proc_id++;
487 :
488 88545 : ids_to_send[proc_id].push_back(dof_id);
489 :
490 : // Note this invalidates the data in nonlocal_pattern
491 88545 : rows_to_send[proc_id].push_back(std::move(row));
492 :
493 : // Might as well remove it since it's invalidated anyway
494 57794 : it = nonlocal_pattern.erase(it);
495 : }
496 :
497 2828 : std::map<processor_id_type, std::vector<dof_id_type>> received_ids_map;
498 :
499 : auto ids_action_functor =
500 : [& received_ids_map]
501 : (processor_id_type pid,
502 1474 : const std::vector<dof_id_type> & received_ids)
503 : {
504 1474 : received_ids_map.emplace(pid, received_ids);
505 5469 : };
506 :
507 4680 : Parallel::push_parallel_vector_data(this->comm(), ids_to_send,
508 : ids_action_functor);
509 :
510 : auto rows_action_functor =
511 685 : [this,
512 : & received_ids_map,
513 : local_first_dof]
514 : (processor_id_type pid,
515 183061 : const std::vector<Row> & received_rows)
516 : {
517 2263 : const std::vector<dof_id_type> & received_ids = libmesh_map_find(received_ids_map, pid);
518 :
519 1578 : std::size_t n_rows = received_rows.size();
520 789 : libmesh_assert_equal_to(n_rows, received_ids.size());
521 :
522 90808 : for (auto i : IntRange<std::size_t>(0, n_rows))
523 : {
524 88545 : const auto r = received_ids[i];
525 30751 : libmesh_assert(dof_map.local_index(r));
526 :
527 88545 : const auto my_r = r - local_first_dof;
528 :
529 61502 : auto & their_row = received_rows[i];
530 :
531 61502 : auto & my_row = sparsity_pattern[my_r];
532 :
533 : // They wouldn't have sent an empty row
534 30751 : libmesh_assert(!their_row.empty());
535 :
536 : // We can end up with an empty row on a dof that touches our
537 : // inactive elements but not our active ones
538 88545 : if (my_row.empty())
539 : {
540 723 : my_row.assign (their_row.begin(), their_row.end());
541 : }
542 : else
543 : {
544 57312 : my_row.insert (my_row.end(),
545 : their_row.begin(),
546 122040 : their_row.end());
547 :
548 : // We cannot use SparsityPattern::sort_row() here because it expects
549 : // the [begin,middle) [middle,end) to be non-overlapping. This is not
550 : // necessarily the case here, so use std::sort()
551 87822 : std::sort (my_row.begin(), my_row.end());
552 :
553 87822 : my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
554 : }
555 :
556 : }
557 4844 : };
558 :
559 4680 : Parallel::push_parallel_vector_data(this->comm(), rows_to_send,
560 : rows_action_functor);
561 :
562 : // We should have sent everything at this point.
563 1414 : libmesh_assert (nonlocal_pattern.empty());
564 :
565 : // assert these are empty because std::vector::resize will only append the specified element value
566 : // if the new size is greater than the current size. Elements whose indices are less than the
567 : // current size are untouched
568 1414 : libmesh_assert(n_nz.empty());
569 1414 : libmesh_assert(n_oz.empty());
570 4680 : n_nz.resize (n_dofs_on_proc, 0);
571 4680 : n_oz.resize (n_dofs_on_proc, 0);
572 :
573 4680 : const dof_id_type first_dof_on_proc = dof_map.first_dof();
574 1414 : const dof_id_type end_dof_on_proc = dof_map.end_dof();
575 :
576 3223826 : for (dof_id_type i=0; i<n_dofs_on_proc; i++)
577 : {
578 : // Get the row of the sparsity pattern
579 2383630 : SparsityPattern::Row & row = sparsity_pattern[i];
580 :
581 126629120 : for (const auto & df : row)
582 123409974 : if ((df < first_dof_on_proc) || (df >= end_dof_on_proc))
583 10542453 : n_oz[i]++;
584 : else
585 145327671 : n_nz[i]++;
586 :
587 864778 : libmesh_assert(n_nz[i] <= n_dofs_on_proc);
588 :
589 : // If we're not building a full sparsity pattern, then we want
590 : // to avoid overcounting these entries as much as possible.
591 3219146 : if (!need_full_sparsity_pattern)
592 864778 : row.clear();
593 : }
594 4680 : }
595 :
596 :
597 0 : void Build::apply_extra_sparsity_object(SparsityPattern::AugmentSparsityPattern & asp)
598 : {
599 0 : asp.augment_sparsity_pattern (sparsity_pattern, n_nz, n_oz);
600 0 : }
601 :
602 :
603 0 : std::size_t Build::n_nonzeros() const
604 : {
605 : // At some point I'll remember that "C++17" compilers don't always
606 : // come with complete C++17 standard libraries.
607 : // std::size_t total_nonzeros = std::reduce(n_nz.begin(), n_nz.end(), std::size_t(0));
608 : // total_nonzeros += std::reduce(n_oz.begin(), n_oz.end(), std::size_t(0));
609 :
610 0 : std::size_t total_nonzeros = 0;
611 0 : for (auto nnzi : n_nz)
612 0 : total_nonzeros += nnzi;
613 0 : for (auto nozi : n_oz)
614 0 : total_nonzeros += nozi;
615 :
616 0 : this->comm().sum(total_nonzeros);
617 0 : return total_nonzeros;
618 : }
619 :
620 :
621 : } // namespace SparsityPattern
622 : } // namespace libMesh
|
