1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_affine_constraints_templates_h
17 #define dealii_affine_constraints_templates_h
18
19 #include <deal.II/base/config.h>
20
21 #include <deal.II/base/cuda_size.h>
22 #include <deal.II/base/memory_consumption.h>
23 #include <deal.II/base/table.h>
24 #include <deal.II/base/thread_local_storage.h>
25
26 #include <deal.II/lac/affine_constraints.h>
27 #include <deal.II/lac/block_sparse_matrix.h>
28 #include <deal.II/lac/block_sparse_matrix_ez.h>
29 #include <deal.II/lac/block_sparsity_pattern.h>
30 #include <deal.II/lac/block_vector.h>
31 #include <deal.II/lac/chunk_sparse_matrix.h>
32 #include <deal.II/lac/diagonal_matrix.h>
33 #include <deal.II/lac/dynamic_sparsity_pattern.h>
34 #include <deal.II/lac/full_matrix.h>
35 #include <deal.II/lac/la_parallel_block_vector.h>
36 #include <deal.II/lac/la_parallel_vector.h>
37 #include <deal.II/lac/la_vector.h>
38 #include <deal.II/lac/matrix_block.h>
39 #include <deal.II/lac/petsc_block_sparse_matrix.h>
40 #include <deal.II/lac/petsc_block_vector.h>
41 #include <deal.II/lac/petsc_sparse_matrix.h>
42 #include <deal.II/lac/petsc_vector.h>
43 #include <deal.II/lac/sparse_matrix.h>
44 #include <deal.II/lac/sparse_matrix_ez.h>
45 #include <deal.II/lac/sparsity_pattern.h>
46 #include <deal.II/lac/trilinos_block_sparse_matrix.h>
47 #include <deal.II/lac/trilinos_parallel_block_vector.h>
48 #include <deal.II/lac/trilinos_sparse_matrix.h>
49 #include <deal.II/lac/trilinos_vector.h>
50
51 #include <boost/serialization/complex.hpp>
52 #include <boost/serialization/utility.hpp>
53
54 #include <algorithm>
55 #include <complex>
56 #include <iomanip>
57 #include <numeric>
58 #include <ostream>
59 #include <set>
60
61 DEAL_II_NAMESPACE_OPEN
62
63
64
65 template <typename number>
66 void
copy_from(const AffineConstraints<number> & other)67 AffineConstraints<number>::copy_from(const AffineConstraints<number> &other)
68 {
69 lines = other.lines;
70 lines_cache = other.lines_cache;
71 local_lines = other.local_lines;
72 sorted = other.sorted;
73 }
74
75
76
77 template <typename number>
78 bool
79 AffineConstraints<number>::ConstraintLine::
80 operator<(const ConstraintLine &a) const
81 {
82 return index < a.index;
83 }
84
85
86
87 template <typename number>
88 bool
89 AffineConstraints<number>::ConstraintLine::
90 operator==(const ConstraintLine &a) const
91 {
92 return index == a.index;
93 }
94
95
96
97 template <typename number>
98 std::size_t
memory_consumption()99 AffineConstraints<number>::ConstraintLine::memory_consumption() const
100 {
101 return (MemoryConsumption::memory_consumption(index) +
102 MemoryConsumption::memory_consumption(entries) +
103 MemoryConsumption::memory_consumption(inhomogeneity));
104 }
105
106
107
108 template <typename number>
109 const typename AffineConstraints<number>::LineRange
get_lines()110 AffineConstraints<number>::get_lines() const
111 {
112 return boost::make_iterator_range(lines.begin(), lines.end());
113 }
114
115
116
117 template <typename number>
118 bool
is_consistent_in_parallel(const std::vector<IndexSet> & locally_owned_dofs,const IndexSet & locally_active_dofs,const MPI_Comm mpi_communicator,const bool verbose)119 AffineConstraints<number>::is_consistent_in_parallel(
120 const std::vector<IndexSet> &locally_owned_dofs,
121 const IndexSet & locally_active_dofs,
122 const MPI_Comm mpi_communicator,
123 const bool verbose) const
124 {
125 // Helper to return a ConstraintLine object that belongs to row @p row.
126 // If @p row is not constrained or not stored locally, return an empty
127 // constraint object that would correspond to a zero constraints
128 auto get_line = [&](const size_type line_n) -> ConstraintLine {
129 const size_type line_index = calculate_line_index(line_n);
130 if (line_index >= lines_cache.size() ||
131 lines_cache[line_index] == numbers::invalid_size_type)
132 {
133 const ConstraintLine empty = {line_n, {}, 0.0};
134 return empty;
135 }
136 else
137 return lines[lines_cache[line_index]];
138 };
139
140 // identify non-owned rows and send to owner:
141 std::map<unsigned int, std::vector<ConstraintLine>> to_send;
142
143 const unsigned int myid =
144 dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
145 const unsigned int nproc =
146 dealii::Utilities::MPI::n_mpi_processes(mpi_communicator);
147
148 // We will send all locally active dofs that are not locally owned for
149 // checking. Note that we allow constraints to differ on locally_relevant (and
150 // not active) DoFs.
151 IndexSet non_owned = locally_active_dofs;
152 non_owned.subtract_set(locally_owned_dofs[myid]);
153 for (unsigned int owner = 0; owner < nproc; ++owner)
154 {
155 // find all lines to send to @p owner
156 IndexSet indices_to_send = non_owned & locally_owned_dofs[owner];
157 for (const auto row_idx : indices_to_send)
158 {
159 to_send[owner].emplace_back(get_line(row_idx));
160 }
161 }
162
163 std::map<unsigned int, std::vector<ConstraintLine>> received =
164 Utilities::MPI::some_to_some(mpi_communicator, to_send);
165
166 unsigned int inconsistent = 0;
167
168 // from each processor:
169 for (const auto &kv : received)
170 {
171 // for each incoming line:
172 for (const auto &lineit : kv.second)
173 {
174 const ConstraintLine reference = get_line(lineit.index);
175
176 if (lineit.inhomogeneity != reference.inhomogeneity)
177 {
178 ++inconsistent;
179
180 if (verbose)
181 std::cout << "Proc " << myid << " got line " << lineit.index
182 << " from " << kv.first << " inhomogeneity "
183 << lineit.inhomogeneity
184 << " != " << reference.inhomogeneity << std::endl;
185 }
186 else if (lineit.entries != reference.entries)
187 {
188 ++inconsistent;
189 if (verbose)
190 std::cout << "Proc " << myid << " got line " << lineit.index
191 << " from " << kv.first << " with wrong values!"
192 << std::endl;
193 }
194 }
195 }
196
197 const unsigned int total =
198 Utilities::MPI::sum(inconsistent, mpi_communicator);
199 if (verbose && total > 0 && myid == 0)
200 std::cout << total << " inconsistent lines discovered!" << std::endl;
201 return total == 0;
202 }
203
204
205
206 template <typename number>
207 void
add_lines(const std::set<size_type> & lines)208 AffineConstraints<number>::add_lines(const std::set<size_type> &lines)
209 {
210 for (const size_type &line : lines)
211 add_line(line);
212 }
213
214
215
216 template <typename number>
217 void
add_lines(const std::vector<bool> & lines)218 AffineConstraints<number>::add_lines(const std::vector<bool> &lines)
219 {
220 for (size_type i = 0; i < lines.size(); ++i)
221 if (lines[i] == true)
222 add_line(i);
223 }
224
225
226
227 template <typename number>
228 void
add_lines(const IndexSet & lines)229 AffineConstraints<number>::add_lines(const IndexSet &lines)
230 {
231 for (size_type i = 0; i < lines.n_elements(); ++i)
232 add_line(lines.nth_index_in_set(i));
233 }
234
235
236
237 template <typename number>
238 void
add_entries(const size_type line_n,const std::vector<std::pair<size_type,number>> & col_val_pairs)239 AffineConstraints<number>::add_entries(
240 const size_type line_n,
241 const std::vector<std::pair<size_type, number>> &col_val_pairs)
242 {
243 Assert(sorted == false, ExcMatrixIsClosed());
244 Assert(is_constrained(line_n), ExcLineInexistant(line_n));
245
246 ConstraintLine &line = lines[lines_cache[calculate_line_index(line_n)]];
247 Assert(line.index == line_n, ExcInternalError());
248
249 // if in debug mode, check whether an entry for this column already
250 // exists and if its the same as the one entered at present
251 //
252 // in any case: skip this entry if an entry for this column already
253 // exists, since we don't want to enter it twice
254 for (const std::pair<size_type, number> &col_val_pair : col_val_pairs)
255 {
256 Assert(line_n != col_val_pair.first,
257 ExcMessage("Can't constrain a degree of freedom to itself"));
258 bool entry_exists = false;
259 for (const std::pair<size_type, number> &entry : line.entries)
260 if (entry.first == col_val_pair.first)
261 {
262 // entry exists, break innermost loop
263 Assert(entry.second == col_val_pair.second,
264 ExcEntryAlreadyExists(line_n,
265 col_val_pair.first,
266 entry.second,
267 col_val_pair.second));
268 entry_exists = true;
269 break;
270 }
271
272 if (entry_exists == false)
273 line.entries.push_back(col_val_pair);
274 }
275 }
276
277
278
279 template <typename number>
280 void
add_selected_constraints(const AffineConstraints & constraints,const IndexSet & filter)281 AffineConstraints<number>::add_selected_constraints(
282 const AffineConstraints &constraints,
283 const IndexSet & filter)
284 {
285 if (constraints.n_constraints() == 0)
286 return;
287
288 Assert(filter.size() > constraints.lines.back().index,
289 ExcMessage(
290 "The filter must be larger than the given constraints object."));
291 for (const ConstraintLine &line : constraints.lines)
292 if (filter.is_element(line.index))
293 {
294 const size_type row = filter.index_within_set(line.index);
295 add_line(row);
296 set_inhomogeneity(row, line.inhomogeneity);
297 for (const std::pair<size_type, number> &entry : line.entries)
298 if (filter.is_element(entry.first))
299 add_entry(row, filter.index_within_set(entry.first), entry.second);
300 }
301 }
302
303
304
305 template <typename number>
306 void
close()307 AffineConstraints<number>::close()
308 {
309 if (sorted == true)
310 return;
311
312 // sort the lines
313 std::sort(lines.begin(), lines.end());
314
315 // update list of pointers and give the vector a sharp size since we
316 // won't modify the size any more after this point.
317 {
318 std::vector<size_type> new_lines(lines_cache.size(),
319 numbers::invalid_size_type);
320 size_type counter = 0;
321 for (const ConstraintLine &line : lines)
322 {
323 new_lines[calculate_line_index(line.index)] = counter;
324 ++counter;
325 }
326 std::swap(lines_cache, new_lines);
327 }
328
329 // in debug mode: check whether we really set the pointers correctly.
330 for (size_type i = 0; i < lines_cache.size(); ++i)
331 if (lines_cache[i] != numbers::invalid_size_type)
332 Assert(i == calculate_line_index(lines[lines_cache[i]].index),
333 ExcInternalError());
334
335 // first, strip zero entries, as we have to do that only once
336 for (ConstraintLine &line : lines)
337 // first remove zero entries. that would mean that in the linear
338 // constraint for a node, x_i = ax_1 + bx_2 + ..., another node times 0
339 // appears. obviously, 0*something can be omitted
340 line.entries.erase(std::remove_if(line.entries.begin(),
341 line.entries.end(),
342 [](
343 const std::pair<size_type, number> &p) {
344 return p.second == number(0.);
345 }),
346 line.entries.end());
347
348
349
350 #ifdef DEBUG
351 // In debug mode we are computing an estimate for the maximum number
352 // of constraints so that we can bail out if there is a cycle in the
353 // constraints (which is easier than searching for cycles in the graph).
354 //
355 // Let us figure out the largest dof index. This is an upper bound for the
356 // number of constraints because it is an approximation for the number of dofs
357 // in our system.
358 size_type largest_idx = 0;
359 for (const ConstraintLine &line : lines)
360 for (const std::pair<size_type, number> &entry : line.entries)
361 largest_idx = std::max(largest_idx, entry.first);
362 #endif
363
364 // replace references to dofs that are themselves constrained. note that
365 // because we may replace references to other dofs that may themselves be
366 // constrained to third ones, we have to iterate over all this until we
367 // replace no chains of constraints any more
368 //
369 // the iteration replaces references to constrained degrees of freedom by
370 // second-order references. for example if x3=x0/2+x2/2 and x2=x0/2+x1/2,
371 // then the new list will be x3=x0/2+x0/4+x1/4. note that x0 appear
372 // twice. we will throw this duplicate out in the following step, where
373 // we sort the list so that throwing out duplicates becomes much more
374 // efficient. also, we have to do it only once, rather than in each
375 // iteration
376 size_type iteration = 0;
377 while (true)
378 {
379 bool chained_constraint_replaced = false;
380
381 for (ConstraintLine &line : lines)
382 {
383 #ifdef DEBUG
384 // we need to keep track of how many replacements we do in this line,
385 // because we can end up in a cycle A->B->C->A without the number of
386 // entries growing.
387 size_type n_replacements = 0;
388 #endif
389
390 // loop over all entries of this line (including ones that we
391 // have appended in this go around) and see whether they are
392 // further constrained. ignore elements that we don't store on
393 // the current processor
394 size_type entry = 0;
395 while (entry < line.entries.size())
396 if (((local_lines.size() == 0) ||
397 (local_lines.is_element(line.entries[entry].first))) &&
398 is_constrained(line.entries[entry].first))
399 {
400 // ok, this entry is further constrained:
401 chained_constraint_replaced = true;
402
403 // look up the chain of constraints for this entry
404 const size_type dof_index = line.entries[entry].first;
405 const number weight = line.entries[entry].second;
406
407 Assert(dof_index != line.index,
408 ExcMessage("Cycle in constraints detected!"));
409
410 const ConstraintLine &constrained_line =
411 lines[lines_cache[calculate_line_index(dof_index)]];
412 Assert(constrained_line.index == dof_index, ExcInternalError());
413
414 // now we have to replace an entry by its expansion. we do
415 // that by overwriting the entry by the first entry of the
416 // expansion and adding the remaining ones to the end,
417 // where we will later process them once more
418 //
419 // we can of course only do that if the DoF that we are
420 // currently handle is constrained by a linear combination
421 // of other dofs:
422 if (constrained_line.entries.size() > 0)
423 {
424 for (size_type i = 0; i < constrained_line.entries.size();
425 ++i)
426 Assert(dof_index != constrained_line.entries[i].first,
427 ExcMessage("Cycle in constraints detected!"));
428
429 // replace first entry, then tack the rest to the end
430 // of the list
431 line.entries[entry] = std::pair<size_type, number>(
432 constrained_line.entries[0].first,
433 constrained_line.entries[0].second * weight);
434
435 for (size_type i = 1; i < constrained_line.entries.size();
436 ++i)
437 line.entries.emplace_back(
438 constrained_line.entries[i].first,
439 constrained_line.entries[i].second * weight);
440
441 #ifdef DEBUG
442 // keep track of how many entries we replace in this
443 // line. If we do more than there are constraints or
444 // dofs in our system, we must have a cycle.
445 ++n_replacements;
446 Assert(n_replacements / 2 < largest_idx,
447 ExcMessage("Cycle in constraints detected!"));
448 if (n_replacements / 2 >= largest_idx)
449 return; // this enables us to test for this Exception.
450 #endif
451 }
452 else
453 // the DoF that we encountered is not constrained by a
454 // linear combination of other dofs but is equal to just
455 // the inhomogeneity (i.e. its chain of entries is
456 // empty). in that case, we can't just overwrite the
457 // current entry, but we have to actually eliminate it
458 {
459 line.entries.erase(line.entries.begin() + entry);
460 }
461
462 line.inhomogeneity += constrained_line.inhomogeneity * weight;
463
464 // now that we're here, do not increase index by one but
465 // rather make another pass for the present entry because
466 // we have replaced the present entry by another one, or
467 // because we have deleted it and shifted all following
468 // ones one forward
469 }
470 else
471 // entry not further constrained. just move ahead by one
472 ++entry;
473 }
474
475 // if we didn't do anything in this round, then quit the loop
476 if (chained_constraint_replaced == false)
477 break;
478
479 // increase iteration count. note that we should not iterate more
480 // times than there are constraints, since this puts a natural upper
481 // bound on the length of constraint chains
482 ++iteration;
483 Assert(iteration <= lines.size(), ExcInternalError());
484 }
485
486 // finally sort the entries and re-scale them if necessary. in this step,
487 // we also throw out duplicates as mentioned above. moreover, as some
488 // entries might have had zero weights, we replace them by a vector with
489 // sharp sizes.
490 for (ConstraintLine &line : lines)
491 {
492 std::sort(line.entries.begin(),
493 line.entries.end(),
494 [](const std::pair<unsigned int, number> &a,
495 const std::pair<unsigned int, number> &b) -> bool {
496 // Let's use lexicogrpahic ordering with std::abs for number
497 // type (it might be complex valued).
498 return (a.first < b.first) ||
499 (a.first == b.first &&
500 std::abs(a.second) < std::abs(b.second));
501 });
502
503 // loop over the now sorted list and see whether any of the entries
504 // references the same dofs more than once in order to find how many
505 // non-duplicate entries we have. This lets us allocate the correct
506 // amount of memory for the constraint entries.
507 size_type duplicates = 0;
508 for (size_type i = 1; i < line.entries.size(); ++i)
509 if (line.entries[i].first == line.entries[i - 1].first)
510 duplicates++;
511
512 if (duplicates > 0 || line.entries.size() < line.entries.capacity())
513 {
514 typename ConstraintLine::Entries new_entries;
515
516 // if we have no duplicates, copy verbatim the entries. this way,
517 // the final size is of the vector is correct.
518 if (duplicates == 0)
519 new_entries = line.entries;
520 else
521 {
522 // otherwise, we need to go through the list and resolve the
523 // duplicates
524 new_entries.reserve(line.entries.size() - duplicates);
525 new_entries.push_back(line.entries[0]);
526 for (size_type j = 1; j < line.entries.size(); ++j)
527 if (line.entries[j].first == line.entries[j - 1].first)
528 {
529 Assert(new_entries.back().first == line.entries[j].first,
530 ExcInternalError());
531 new_entries.back().second += line.entries[j].second;
532 }
533 else
534 new_entries.push_back(line.entries[j]);
535
536 Assert(new_entries.size() == line.entries.size() - duplicates,
537 ExcInternalError());
538
539 // make sure there are really no duplicates left and that the
540 // list is still sorted
541 for (size_type j = 1; j < new_entries.size(); ++j)
542 {
543 Assert(new_entries[j].first != new_entries[j - 1].first,
544 ExcInternalError());
545 Assert(new_entries[j].first > new_entries[j - 1].first,
546 ExcInternalError());
547 }
548 }
549
550 // replace old list of constraints for this dof by the new one
551 line.entries.swap(new_entries);
552 }
553
554 // Finally do the following check: if the sum of weights for the
555 // constraints is close to one, but not exactly one, then rescale all
556 // the weights so that they sum up to 1. this adds a little numerical
557 // stability and avoids all sorts of problems where the actual value
558 // is close to, but not quite what we expected
559 //
560 // the case where the weights don't quite sum up happens when we
561 // compute the interpolation weights "on the fly", i.e. not from
562 // precomputed tables. in this case, the interpolation weights are
563 // also subject to round-off
564 number sum = 0.;
565 for (const std::pair<size_type, number> &entry : line.entries)
566 sum += entry.second;
567 if (std::abs(sum - number(1.)) < 1.e-13)
568 {
569 for (std::pair<size_type, number> &entry : line.entries)
570 entry.second /= sum;
571 line.inhomogeneity /= sum;
572 }
573 } // end of loop over all constraint lines
574
575 #ifdef DEBUG
576 // if in debug mode: check that no dof is constrained to another dof that
577 // is also constrained. exclude dofs from this check whose constraint
578 // lines are not stored on the local processor
579 for (const ConstraintLine &line : lines)
580 for (const std::pair<size_type, number> &entry : line.entries)
581 if ((local_lines.size() == 0) || (local_lines.is_element(entry.first)))
582 {
583 // make sure that entry->first is not the index of a line itself
584 const bool is_circle = is_constrained(entry.first);
585 Assert(is_circle == false,
586 ExcDoFConstrainedToConstrainedDoF(line.index, entry.first));
587 }
588 #endif
589
590 sorted = true;
591 }
592
593
594
595 template <typename number>
596 void
merge(const AffineConstraints<number> & other_constraints,const MergeConflictBehavior merge_conflict_behavior,const bool allow_different_local_lines)597 AffineConstraints<number>::merge(
598 const AffineConstraints<number> &other_constraints,
599 const MergeConflictBehavior merge_conflict_behavior,
600 const bool allow_different_local_lines)
601 {
602 (void)allow_different_local_lines;
603 Assert(allow_different_local_lines ||
604 local_lines == other_constraints.local_lines,
605 ExcMessage(
606 "local_lines for this and the other objects are not the same "
607 "although allow_different_local_lines is false."));
608
609 // store the previous state with respect to sorting
610 const bool object_was_sorted = sorted;
611 sorted = false;
612
613 // first action is to fold into the present object possible constraints
614 // in the second object. we don't strictly need to do this any more since
615 // the AffineConstraints container has learned to deal with chains of
616 // constraints in the close() function, but we have traditionally done
617 // this and it's not overly hard to do.
618 //
619 // for this, loop over all constraints and replace the constraint lines
620 // with a new one where constraints are replaced if necessary.
621 typename ConstraintLine::Entries tmp;
622 for (ConstraintLine &line : lines)
623 {
624 tmp.clear();
625 for (const std::pair<size_type, number> &entry : line.entries)
626 {
627 // if the present dof is not stored, or not constrained, or if we
628 // won't take the constraint from the other object, then simply copy
629 // it over
630 if ((other_constraints.local_lines.size() != 0. &&
631 other_constraints.local_lines.is_element(entry.first) ==
632 false) ||
633 other_constraints.is_constrained(entry.first) == false ||
634 ((merge_conflict_behavior != right_object_wins) &&
635 other_constraints.is_constrained(entry.first) &&
636 this->is_constrained(entry.first)))
637 tmp.push_back(entry);
638 else
639 // otherwise resolve further constraints by replacing the old
640 // entry by a sequence of new entries taken from the other
641 // object, but with multiplied weights
642 {
643 const typename ConstraintLine::Entries *other_entries =
644 other_constraints.get_constraint_entries(entry.first);
645 Assert(other_entries != nullptr, ExcInternalError());
646
647 const number weight = entry.second;
648
649 for (const std::pair<size_type, number> &other_entry :
650 *other_entries)
651 tmp.emplace_back(other_entry.first,
652 other_entry.second * weight);
653
654 line.inhomogeneity +=
655 other_constraints.get_inhomogeneity(entry.first) * weight;
656 }
657 }
658 // finally exchange old and newly resolved line
659 line.entries.swap(tmp);
660 }
661
662 if (local_lines.size() != 0)
663 local_lines.add_indices(other_constraints.local_lines);
664
665 {
666 // do not bother to resize the lines cache exactly since it is pretty
667 // cheap to adjust it along the way.
668 std::fill(lines_cache.begin(),
669 lines_cache.end(),
670 numbers::invalid_size_type);
671
672 // reset lines_cache for our own constraints
673 size_type index = 0;
674 for (const ConstraintLine &line : lines)
675 {
676 const size_type local_line_no = calculate_line_index(line.index);
677 if (local_line_no >= lines_cache.size())
678 lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
679 lines_cache[local_line_no] = index++;
680 }
681
682 // Add other_constraints to lines cache and our list of constraints
683 for (const ConstraintLine &line : other_constraints.lines)
684 {
685 const size_type local_line_no = calculate_line_index(line.index);
686 if (local_line_no >= lines_cache.size())
687 {
688 lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
689 lines.push_back(line);
690 lines_cache[local_line_no] = index++;
691 }
692 else if (lines_cache[local_line_no] == numbers::invalid_size_type)
693 {
694 // there are no constraints for that line yet
695 lines.push_back(line);
696 AssertIndexRange(local_line_no, lines_cache.size());
697 lines_cache[local_line_no] = index++;
698 }
699 else
700 {
701 // we already store that line
702 switch (merge_conflict_behavior)
703 {
704 case no_conflicts_allowed:
705 AssertThrow(false,
706 ExcDoFIsConstrainedFromBothObjects(line.index));
707 break;
708
709 case left_object_wins:
710 // ignore this constraint
711 break;
712
713 case right_object_wins:
714 AssertIndexRange(local_line_no, lines_cache.size());
715 lines[lines_cache[local_line_no]] = line;
716 break;
717
718 default:
719 Assert(false, ExcNotImplemented());
720 }
721 }
722 }
723
724 // check that we set the pointers correctly
725 for (size_type i = 0; i < lines_cache.size(); ++i)
726 if (lines_cache[i] != numbers::invalid_size_type)
727 Assert(i == calculate_line_index(lines[lines_cache[i]].index),
728 ExcInternalError());
729 }
730
731 // if the object was sorted before, then make sure it is so afterward as
732 // well. otherwise leave everything in the unsorted state
733 if (object_was_sorted == true)
734 close();
735 }
736
737
738
739 template <typename number>
740 void
shift(const size_type offset)741 AffineConstraints<number>::shift(const size_type offset)
742 {
743 if (local_lines.size() == 0)
744 lines_cache.insert(lines_cache.begin(), offset, numbers::invalid_size_type);
745 else
746 {
747 // shift local_lines
748 IndexSet new_local_lines(local_lines.size());
749 new_local_lines.add_indices(local_lines, offset);
750 std::swap(local_lines, new_local_lines);
751 }
752
753 for (ConstraintLine &line : lines)
754 {
755 line.index += offset;
756 for (std::pair<size_type, number> &entry : line.entries)
757 entry.first += offset;
758 }
759
760 #ifdef DEBUG
761 // make sure that lines, lines_cache and local_lines
762 // are still linked correctly
763 for (size_type index = 0; index < lines_cache.size(); ++index)
764 Assert(lines_cache[index] == numbers::invalid_size_type ||
765 calculate_line_index(lines[lines_cache[index]].index) == index,
766 ExcInternalError());
767 #endif
768 }
769
770
771
772 template <typename number>
773 void
clear()774 AffineConstraints<number>::clear()
775 {
776 {
777 std::vector<ConstraintLine> tmp;
778 lines.swap(tmp);
779 }
780
781 {
782 std::vector<size_type> tmp;
783 lines_cache.swap(tmp);
784 }
785
786 sorted = false;
787 }
788
789
790
791 template <typename number>
792 void
reinit(const IndexSet & local_constraints)793 AffineConstraints<number>::reinit(const IndexSet &local_constraints)
794 {
795 local_lines = local_constraints;
796
797 // make sure the IndexSet is compressed. Otherwise this can lead to crashes
798 // that are hard to find (only happen in release mode).
799 // see tests/mpi/affine_constraints_crash_01
800 local_lines.compress();
801
802 clear();
803 }
804
805
806
807 template <typename number>
808 bool
is_identity_constrained(const size_type line_n)809 AffineConstraints<number>::is_identity_constrained(const size_type line_n) const
810 {
811 if (is_constrained(line_n) == false)
812 return false;
813
814 const ConstraintLine &line = lines[lines_cache[calculate_line_index(line_n)]];
815 Assert(line.index == line_n, ExcInternalError());
816
817 // return if an entry for this line was found and if it has only one
818 // entry equal to 1.0
819 return (line.entries.size() == 1) && (line.entries[0].second == number(1.0));
820 }
821
822
823 template <typename number>
824 bool
are_identity_constrained(const size_type line_n_1,const size_type line_n_2)825 AffineConstraints<number>::are_identity_constrained(
826 const size_type line_n_1,
827 const size_type line_n_2) const
828 {
829 if (is_constrained(line_n_1) == true)
830 {
831 const ConstraintLine &line =
832 lines[lines_cache[calculate_line_index(line_n_1)]];
833 Assert(line.index == line_n_1, ExcInternalError());
834
835 // return if an entry for this line was found and if it has only one
836 // entry equal to 1.0 and that one is index2
837 return ((line.entries.size() == 1) &&
838 (line.entries[0].first == line_n_2) &&
839 (line.entries[0].second == number(1.0)));
840 }
841 else if (is_constrained(line_n_2) == true)
842 {
843 const ConstraintLine &line =
844 lines[lines_cache[calculate_line_index(line_n_2)]];
845 Assert(line.index == line_n_2, ExcInternalError());
846
847 // return if an entry for this line was found and if it has only one
848 // entry equal to 1.0 and that one is line_n_1
849 return ((line.entries.size() == 1) &&
850 (line.entries[0].first == line_n_1) &&
851 (line.entries[0].second == number(1.0)));
852 }
853 else
854 return false;
855 }
856
857
858
859 template <typename number>
860 typename AffineConstraints<number>::size_type
max_constraint_indirections()861 AffineConstraints<number>::max_constraint_indirections() const
862 {
863 size_type return_value = 0;
864 for (const ConstraintLine &line : lines)
865 // use static cast, since typeof(size)==std::size_t, which is !=
866 // size_type on AIX
867 return_value =
868 std::max(return_value, static_cast<size_type>(line.entries.size()));
869
870 return return_value;
871 }
872
873
874
875 template <typename number>
876 bool
has_inhomogeneities()877 AffineConstraints<number>::has_inhomogeneities() const
878 {
879 for (const ConstraintLine &line : lines)
880 if (line.inhomogeneity != number(0.))
881 return true;
882
883 return false;
884 }
885
886
887
888 template <typename number>
889 void
print(std::ostream & out)890 AffineConstraints<number>::print(std::ostream &out) const
891 {
892 for (const ConstraintLine &line : lines)
893 {
894 // output the list of constraints as pairs of dofs and their weights
895 if (line.entries.size() > 0)
896 {
897 for (const std::pair<size_type, number> &entry : line.entries)
898 out << " " << line.index << " " << entry.first << ": "
899 << entry.second << "\n";
900
901 // print out inhomogeneity.
902 if (line.inhomogeneity != number(0.))
903 out << " " << line.index << ": " << line.inhomogeneity << "\n";
904 }
905 else
906 // but also output something if the constraint simply reads
907 // x[13]=0, i.e. where the right hand side is not a linear
908 // combination of other dofs
909 {
910 if (line.inhomogeneity != number(0.))
911 out << " " << line.index << " = " << line.inhomogeneity << "\n";
912 else
913 out << " " << line.index << " = 0\n";
914 }
915 }
916
917 AssertThrow(out, ExcIO());
918 }
919
920
921
922 template <typename number>
923 void
write_dot(std::ostream & out)924 AffineConstraints<number>::write_dot(std::ostream &out) const
925 {
926 out << "digraph constraints {" << std::endl;
927 for (size_type i = 0; i != lines.size(); ++i)
928 {
929 // same concept as in the previous function
930 if (lines[i].entries.size() > 0)
931 for (size_type j = 0; j < lines[i].entries.size(); ++j)
932 out << " " << lines[i].index << "->" << lines[i].entries[j].first
933 << "; // weight: " << lines[i].entries[j].second << "\n";
934 else
935 out << " " << lines[i].index << "\n";
936 }
937 out << "}" << std::endl;
938 }
939
940
941
942 template <typename number>
943 std::size_t
memory_consumption()944 AffineConstraints<number>::memory_consumption() const
945 {
946 return (MemoryConsumption::memory_consumption(lines) +
947 MemoryConsumption::memory_consumption(lines_cache) +
948 MemoryConsumption::memory_consumption(sorted) +
949 MemoryConsumption::memory_consumption(local_lines));
950 }
951
952
953
954 template <typename number>
955 void
resolve_indices(std::vector<types::global_dof_index> & indices)956 AffineConstraints<number>::resolve_indices(
957 std::vector<types::global_dof_index> &indices) const
958 {
959 const unsigned int indices_size = indices.size();
960 const std::vector<std::pair<types::global_dof_index, number>> *line_ptr;
961 for (unsigned int i = 0; i < indices_size; ++i)
962 {
963 line_ptr = get_constraint_entries(indices[i]);
964 // if the index is constraint, the constraints indices are added to the
965 // indices vector
966 if (line_ptr != nullptr)
967 {
968 const unsigned int line_size = line_ptr->size();
969 for (unsigned int j = 0; j < line_size; ++j)
970 indices.push_back((*line_ptr)[j].first);
971 }
972 }
973
974 // keep only the unique elements
975 std::sort(indices.begin(), indices.end());
976 std::vector<types::global_dof_index>::iterator it;
977 it = std::unique(indices.begin(), indices.end());
978 indices.resize(it - indices.begin());
979 }
980
981
982
983 template <typename number>
984 void
condense(SparsityPattern & sparsity)985 AffineConstraints<number>::condense(SparsityPattern &sparsity) const
986 {
987 Assert(sorted == true, ExcMatrixNotClosed());
988 Assert(sparsity.is_compressed() == false, ExcMatrixIsClosed());
989 Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
990
991 // store for each index whether it must be distributed or not. If entry
992 // is numbers::invalid_unsigned_int, no distribution is necessary.
993 // otherwise, the number states which line in the AffineConstraints object
994 // handles this index
995 std::vector<size_type> distribute(sparsity.n_rows(),
996 numbers::invalid_size_type);
997
998 for (size_type c = 0; c < lines.size(); ++c)
999 distribute[lines[c].index] = c;
1000
1001 const size_type n_rows = sparsity.n_rows();
1002 for (size_type row = 0; row < n_rows; ++row)
1003 {
1004 if (distribute[row] == numbers::invalid_size_type)
1005 {
1006 // regular line. loop over cols all valid cols. note that this
1007 // changes the line we are presently working on: we add additional
1008 // entries. these are put to the end of the row. however, as
1009 // constrained nodes cannot be constrained to other constrained
1010 // nodes, nothing will happen if we run into these added nodes, as
1011 // they can't be distributed further. we might store the position of
1012 // the last old entry and stop work there, but since operating on
1013 // the newly added ones only takes two comparisons (column index
1014 // valid, distribute[column] necessarily
1015 // ==numbers::invalid_size_type), it is cheaper to not do so and
1016 // run right until the end of the line
1017 for (SparsityPattern::iterator entry = sparsity.begin(row);
1018 ((entry != sparsity.end(row)) && entry->is_valid_entry());
1019 ++entry)
1020 {
1021 const size_type column = entry->column();
1022
1023 if (distribute[column] != numbers::invalid_size_type)
1024 {
1025 // distribute entry at regular row @p{row} and irregular
1026 // column sparsity.colnums[j]
1027 for (size_type q = 0;
1028 q != lines[distribute[column]].entries.size();
1029 ++q)
1030 sparsity.add(row,
1031 lines[distribute[column]].entries[q].first);
1032 }
1033 }
1034 }
1035 else
1036 // row must be distributed. note that here the present row is not
1037 // touched (unlike above)
1038 {
1039 for (SparsityPattern::iterator entry = sparsity.begin(row);
1040 (entry != sparsity.end(row)) && entry->is_valid_entry();
1041 ++entry)
1042 {
1043 const size_type column = entry->column();
1044 if (distribute[column] == numbers::invalid_size_type)
1045 // distribute entry at irregular row @p{row} and regular
1046 // column sparsity.colnums[j]
1047 for (size_type q = 0;
1048 q != lines[distribute[row]].entries.size();
1049 ++q)
1050 sparsity.add(lines[distribute[row]].entries[q].first, column);
1051 else
1052 // distribute entry at irregular row @p{row} and irregular
1053 // column sparsity.get_column_numbers()[j]
1054 for (size_type p = 0;
1055 p != lines[distribute[row]].entries.size();
1056 ++p)
1057 for (size_type q = 0;
1058 q != lines[distribute[column]].entries.size();
1059 ++q)
1060 sparsity.add(lines[distribute[row]].entries[p].first,
1061 lines[distribute[column]].entries[q].first);
1062 }
1063 }
1064 }
1065
1066 sparsity.compress();
1067 }
1068
1069
1070
1071 template <typename number>
1072 void
condense(BlockSparsityPattern & sparsity)1073 AffineConstraints<number>::condense(BlockSparsityPattern &sparsity) const
1074 {
1075 Assert(sorted == true, ExcMatrixNotClosed());
1076 Assert(sparsity.is_compressed() == false, ExcMatrixIsClosed());
1077 Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
1078 Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
1079 Assert(sparsity.get_column_indices() == sparsity.get_row_indices(),
1080 ExcNotQuadratic());
1081
1082 const BlockIndices &index_mapping = sparsity.get_column_indices();
1083
1084 const size_type n_blocks = sparsity.n_block_rows();
1085
1086 // store for each index whether it must be distributed or not. If entry
1087 // is numbers::invalid_unsigned_int, no distribution is necessary.
1088 // otherwise, the number states which line in the AffineConstraints object
1089 // handles this index
1090 std::vector<size_type> distribute(sparsity.n_rows(),
1091 numbers::invalid_size_type);
1092
1093 for (size_type c = 0; c < lines.size(); ++c)
1094 distribute[lines[c].index] = c;
1095
1096 const size_type n_rows = sparsity.n_rows();
1097 for (size_type row = 0; row < n_rows; ++row)
1098 {
1099 // get index of this row within the blocks
1100 const std::pair<size_type, size_type> block_index =
1101 index_mapping.global_to_local(row);
1102 const size_type block_row = block_index.first;
1103
1104 if (distribute[row] == numbers::invalid_size_type)
1105 // regular line. loop over all columns and see whether this column
1106 // must be distributed
1107 {
1108 // to loop over all entries in this row, we have to loop over all
1109 // blocks in this blockrow and the corresponding row therein
1110 for (size_type block_col = 0; block_col < n_blocks; ++block_col)
1111 {
1112 const SparsityPattern &block_sparsity =
1113 sparsity.block(block_row, block_col);
1114
1115 for (SparsityPattern::const_iterator entry =
1116 block_sparsity.begin(block_index.second);
1117 (entry != block_sparsity.end(block_index.second)) &&
1118 entry->is_valid_entry();
1119 ++entry)
1120 {
1121 const size_type global_col =
1122 index_mapping.local_to_global(block_col, entry->column());
1123
1124 if (distribute[global_col] != numbers::invalid_size_type)
1125 // distribute entry at regular row @p{row} and
1126 // irregular column global_col
1127 {
1128 for (size_type q = 0;
1129 q != lines[distribute[global_col]].entries.size();
1130 ++q)
1131 sparsity.add(
1132 row, lines[distribute[global_col]].entries[q].first);
1133 }
1134 }
1135 }
1136 }
1137 else
1138 {
1139 // row must be distributed. split the whole row into the chunks
1140 // defined by the blocks
1141 for (size_type block_col = 0; block_col < n_blocks; ++block_col)
1142 {
1143 const SparsityPattern &block_sparsity =
1144 sparsity.block(block_row, block_col);
1145
1146 for (SparsityPattern::const_iterator entry =
1147 block_sparsity.begin(block_index.second);
1148 (entry != block_sparsity.end(block_index.second)) &&
1149 entry->is_valid_entry();
1150 ++entry)
1151 {
1152 const size_type global_col =
1153 index_mapping.local_to_global(block_col, entry->column());
1154
1155 if (distribute[global_col] == numbers::invalid_size_type)
1156 // distribute entry at irregular row @p{row} and
1157 // regular column global_col.
1158 {
1159 for (size_type q = 0;
1160 q != lines[distribute[row]].entries.size();
1161 ++q)
1162 sparsity.add(lines[distribute[row]].entries[q].first,
1163 global_col);
1164 }
1165 else
1166 // distribute entry at irregular row @p{row} and
1167 // irregular column @p{global_col}
1168 {
1169 for (size_type p = 0;
1170 p != lines[distribute[row]].entries.size();
1171 ++p)
1172 for (size_type q = 0;
1173 q != lines[distribute[global_col]].entries.size();
1174 ++q)
1175 sparsity.add(
1176 lines[distribute[row]].entries[p].first,
1177 lines[distribute[global_col]].entries[q].first);
1178 }
1179 }
1180 }
1181 }
1182 }
1183
1184 sparsity.compress();
1185 }
1186
1187
1188
1189 template <typename number>
1190 void
condense(DynamicSparsityPattern & sparsity)1191 AffineConstraints<number>::condense(DynamicSparsityPattern &sparsity) const
1192 {
1193 Assert(sorted == true, ExcMatrixNotClosed());
1194 Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
1195
1196 // store for each index whether it must be distributed or not. If entry
1197 // is numbers::invalid_unsigned_int, no distribution is necessary.
1198 // otherwise, the number states which line in the constraints object
1199 // handles this index
1200 std::vector<size_type> distribute(sparsity.n_rows(),
1201 numbers::invalid_size_type);
1202
1203 for (size_type c = 0; c < lines.size(); ++c)
1204 distribute[lines[c].index] = c;
1205
1206 const size_type n_rows = sparsity.n_rows();
1207 for (size_type row = 0; row < n_rows; ++row)
1208 {
1209 if (distribute[row] == numbers::invalid_size_type)
1210 // regular line. loop over cols. note that as we proceed to
1211 // distribute cols, the loop may get longer
1212 for (size_type j = 0; j < sparsity.row_length(row); ++j)
1213 {
1214 const size_type column = sparsity.column_number(row, j);
1215
1216 if (distribute[column] != numbers::invalid_size_type)
1217 {
1218 // distribute entry at regular row @p{row} and irregular
1219 // column column. note that this changes the line we are
1220 // presently working on: we add additional entries. if we
1221 // add another entry at a column behind the present one, we
1222 // will encounter it later on (but since it can't be
1223 // further constrained, won't have to do anything about
1224 // it). if we add it up front of the present column, we
1225 // will find the present column later on again as it was
1226 // shifted back (again nothing happens, in particular no
1227 // endless loop, as when we encounter it the second time we
1228 // won't be able to add more entries as they all already
1229 // exist, but we do the same work more often than
1230 // necessary, and the loop gets longer), so move the cursor
1231 // one to the right in the case that we add an entry up
1232 // front that did not exist before. check whether it
1233 // existed before by tracking the length of this row
1234 size_type old_rowlength = sparsity.row_length(row);
1235 for (size_type q = 0;
1236 q != lines[distribute[column]].entries.size();
1237 ++q)
1238 {
1239 const size_type new_col =
1240 lines[distribute[column]].entries[q].first;
1241
1242 sparsity.add(row, new_col);
1243
1244 const size_type new_rowlength = sparsity.row_length(row);
1245 if ((new_col < column) && (old_rowlength != new_rowlength))
1246 ++j;
1247 old_rowlength = new_rowlength;
1248 }
1249 }
1250 }
1251 else
1252 // row must be distributed
1253 for (size_type j = 0; j < sparsity.row_length(row); ++j)
1254 {
1255 const size_type column = sparsity.column_number(row, j);
1256
1257 if (distribute[column] == numbers::invalid_size_type)
1258 // distribute entry at irregular row @p{row} and regular
1259 // column sparsity.colnums[j]
1260 for (size_type q = 0; q != lines[distribute[row]].entries.size();
1261 ++q)
1262 sparsity.add(lines[distribute[row]].entries[q].first, column);
1263 else
1264 // distribute entry at irregular row @p{row} and irregular
1265 // column sparsity.get_column_numbers()[j]
1266 for (size_type p = 0; p != lines[distribute[row]].entries.size();
1267 ++p)
1268 for (size_type q = 0;
1269 q != lines[distribute[sparsity.column_number(row, j)]]
1270 .entries.size();
1271 ++q)
1272 sparsity.add(lines[distribute[row]].entries[p].first,
1273 lines[distribute[sparsity.column_number(row, j)]]
1274 .entries[q]
1275 .first);
1276 }
1277 }
1278 }
1279
1280
1281
1282 template <typename number>
1283 void
condense(BlockDynamicSparsityPattern & sparsity)1284 AffineConstraints<number>::condense(BlockDynamicSparsityPattern &sparsity) const
1285 {
1286 Assert(sorted == true, ExcMatrixNotClosed());
1287 Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
1288 Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
1289 Assert(sparsity.get_column_indices() == sparsity.get_row_indices(),
1290 ExcNotQuadratic());
1291
1292 const BlockIndices &index_mapping = sparsity.get_column_indices();
1293
1294 const size_type n_blocks = sparsity.n_block_rows();
1295
1296 // store for each index whether it must be distributed or not. If entry
1297 // is numbers::invalid_unsigned_int, no distribution is necessary.
1298 // otherwise, the number states which line in the constraints object
1299 // handles this index
1300 std::vector<size_type> distribute(sparsity.n_rows(),
1301 numbers::invalid_size_type);
1302
1303 for (size_type c = 0; c < lines.size(); ++c)
1304 distribute[lines[c].index] = static_cast<signed int>(c);
1305
1306 const size_type n_rows = sparsity.n_rows();
1307 for (size_type row = 0; row < n_rows; ++row)
1308 {
1309 // get index of this row within the blocks
1310 const std::pair<size_type, size_type> block_index =
1311 index_mapping.global_to_local(row);
1312 const size_type block_row = block_index.first;
1313 const size_type local_row = block_index.second;
1314
1315 if (distribute[row] == numbers::invalid_size_type)
1316 // regular line. loop over all columns and see whether this column
1317 // must be distributed. note that as we proceed to distribute cols,
1318 // the loop over cols may get longer.
1319 //
1320 // don't try to be clever here as in the algorithm for the
1321 // DynamicSparsityPattern, as that would be much more
1322 // complicated here. after all, we know that compressed patterns
1323 // are inefficient...
1324 {
1325 // to loop over all entries in this row, we have to loop over all
1326 // blocks in this blockrow and the corresponding row therein
1327 for (size_type block_col = 0; block_col < n_blocks; ++block_col)
1328 {
1329 const DynamicSparsityPattern &block_sparsity =
1330 sparsity.block(block_row, block_col);
1331
1332 for (size_type j = 0; j < block_sparsity.row_length(local_row);
1333 ++j)
1334 {
1335 const size_type global_col = index_mapping.local_to_global(
1336 block_col, block_sparsity.column_number(local_row, j));
1337
1338 if (distribute[global_col] != numbers::invalid_size_type)
1339 // distribute entry at regular row @p{row} and
1340 // irregular column global_col
1341 {
1342 for (size_type q = 0;
1343 q != lines[distribute[global_col]].entries.size();
1344 ++q)
1345 sparsity.add(
1346 row, lines[distribute[global_col]].entries[q].first);
1347 }
1348 }
1349 }
1350 }
1351 else
1352 {
1353 // row must be distributed. split the whole row into the chunks
1354 // defined by the blocks
1355 for (size_type block_col = 0; block_col < n_blocks; ++block_col)
1356 {
1357 const DynamicSparsityPattern &block_sparsity =
1358 sparsity.block(block_row, block_col);
1359
1360 for (size_type j = 0; j < block_sparsity.row_length(local_row);
1361 ++j)
1362 {
1363 const size_type global_col = index_mapping.local_to_global(
1364 block_col, block_sparsity.column_number(local_row, j));
1365
1366 if (distribute[global_col] == numbers::invalid_size_type)
1367 // distribute entry at irregular row @p{row} and
1368 // regular column global_col.
1369 {
1370 for (size_type q = 0;
1371 q != lines[distribute[row]].entries.size();
1372 ++q)
1373 sparsity.add(lines[distribute[row]].entries[q].first,
1374 global_col);
1375 }
1376 else
1377 // distribute entry at irregular row @p{row} and
1378 // irregular column @p{global_col}
1379 {
1380 for (size_type p = 0;
1381 p != lines[distribute[row]].entries.size();
1382 ++p)
1383 for (size_type q = 0;
1384 q != lines[distribute[global_col]].entries.size();
1385 ++q)
1386 sparsity.add(
1387 lines[distribute[row]].entries[p].first,
1388 lines[distribute[global_col]].entries[q].first);
1389 }
1390 }
1391 }
1392 }
1393 }
1394 }
1395
1396
1397
1398 template <typename number>
1399 void
condense(SparseMatrix<number> & uncondensed)1400 AffineConstraints<number>::condense(SparseMatrix<number> &uncondensed) const
1401 {
1402 Vector<number> dummy(0);
1403 condense(uncondensed, dummy);
1404 }
1405
1406
1407
1408 template <typename number>
1409 void
condense(BlockSparseMatrix<number> & uncondensed)1410 AffineConstraints<number>::condense(
1411 BlockSparseMatrix<number> &uncondensed) const
1412 {
1413 BlockVector<number> dummy(0);
1414 condense(uncondensed, dummy);
1415 }
1416
1417
1418
1419 template <typename number>
1420 template <class VectorType>
1421 void
condense(const VectorType & vec_ghosted,VectorType & vec)1422 AffineConstraints<number>::condense(const VectorType &vec_ghosted,
1423 VectorType & vec) const
1424 {
1425 Assert(sorted == true, ExcMatrixNotClosed());
1426
1427 // if this is called with different arguments, we need to copy the data
1428 // over:
1429 if (&vec != &vec_ghosted)
1430 vec = vec_ghosted;
1431
1432 // distribute all entries, and set them to zero. do so in two loops
1433 // because in the first one we need to add to elements and in the second
1434 // one we need to set elements to zero. for parallel vectors, this can
1435 // only work if we can put a compress() in between, but we don't want to
1436 // call compress() twice per entry
1437 for (const ConstraintLine &line : lines)
1438 {
1439 // in case the constraint is inhomogeneous, this function is not
1440 // appropriate. Throw an exception.
1441 Assert(line.inhomogeneity == number(0.),
1442 ExcMessage("Inhomogeneous constraint cannot be condensed "
1443 "without any matrix specified."));
1444
1445 const typename VectorType::value_type old_value = vec_ghosted(line.index);
1446 for (const std::pair<size_type, number> &entry : line.entries)
1447 if (vec.in_local_range(entry.first) == true)
1448 vec(entry.first) +=
1449 (static_cast<typename VectorType::value_type>(old_value) *
1450 entry.second);
1451 }
1452
1453 vec.compress(VectorOperation::add);
1454
1455 for (const ConstraintLine &line : lines)
1456 if (vec.in_local_range(line.index) == true)
1457 vec(line.index) = 0.;
1458
1459 vec.compress(VectorOperation::insert);
1460 }
1461
1462
1463
1464 template <typename number>
1465 template <class VectorType>
1466 void
condense(VectorType & vec)1467 AffineConstraints<number>::condense(VectorType &vec) const
1468 {
1469 condense(vec, vec);
1470 }
1471
1472
1473
1474 template <typename number>
1475 template <class VectorType>
1476 void
condense(SparseMatrix<number> & uncondensed,VectorType & vec)1477 AffineConstraints<number>::condense(SparseMatrix<number> &uncondensed,
1478 VectorType & vec) const
1479 {
1480 // check whether we work on real vectors or we just used a dummy when
1481 // calling the other function above.
1482 const bool use_vectors = vec.size() == 0 ? false : true;
1483
1484 const SparsityPattern &sparsity = uncondensed.get_sparsity_pattern();
1485
1486 Assert(sorted == true, ExcMatrixNotClosed());
1487 Assert(sparsity.is_compressed() == true, ExcMatrixNotClosed());
1488 Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
1489 if (use_vectors == true)
1490 AssertDimension(vec.size(), sparsity.n_rows());
1491
1492 number average_diagonal = 0.;
1493 for (size_type i = 0; i < uncondensed.m(); ++i)
1494 average_diagonal += std::abs(uncondensed.diag_element(i));
1495 average_diagonal /= uncondensed.m();
1496
1497 // store for each index whether it must be distributed or not. If entry
1498 // is invalid_size_type, no distribution is necessary. otherwise, the
1499 // number states which line in the constraints object handles this index
1500 std::vector<size_type> distribute(sparsity.n_rows(),
1501 numbers::invalid_size_type);
1502
1503 for (size_type c = 0; c < lines.size(); ++c)
1504 distribute[lines[c].index] = c;
1505
1506 const size_type n_rows = sparsity.n_rows();
1507 for (size_type row = 0; row < n_rows; ++row)
1508 {
1509 if (distribute[row] == numbers::invalid_size_type)
1510 // regular line. loop over cols
1511 {
1512 for (typename SparseMatrix<number>::iterator entry =
1513 uncondensed.begin(row);
1514 entry != uncondensed.end(row);
1515 ++entry)
1516 {
1517 const size_type column = entry->column();
1518
1519 // end of row reached? this should not happen, since we only
1520 // operate on compressed matrices!
1521 Assert(column != SparsityPattern::invalid_entry,
1522 ExcMatrixNotClosed());
1523
1524 if (distribute[column] != numbers::invalid_size_type)
1525 // distribute entry at regular row @p row and irregular
1526 // column sparsity.get_column_numbers()[j]; set old entry
1527 // to zero
1528 {
1529 for (size_type q = 0;
1530 q != lines[distribute[column]].entries.size();
1531 ++q)
1532 {
1533 // need a temporary variable to avoid errors like no
1534 // known conversion from 'complex<typename
1535 // ProductType<float, double>::type>' to 'const
1536 // complex<float>' for 3rd argument
1537 number v = static_cast<number>(entry->value());
1538 v *= lines[distribute[column]].entries[q].second;
1539 uncondensed.add(
1540 row, lines[distribute[column]].entries[q].first, v);
1541 }
1542
1543 // need to subtract this element from the vector. this
1544 // corresponds to an explicit elimination in the
1545 // respective row of the inhomogeneous constraint in the
1546 // matrix with Gauss elimination
1547 if (use_vectors == true)
1548 vec(row) -= static_cast<number>(entry->value()) *
1549 lines[distribute[column]].inhomogeneity;
1550
1551 // set old value to zero
1552 entry->value() = 0.;
1553 }
1554 }
1555 }
1556 else
1557 // row must be distributed
1558 {
1559 for (typename SparseMatrix<number>::iterator entry =
1560 uncondensed.begin(row);
1561 entry != uncondensed.end(row);
1562 ++entry)
1563 {
1564 const size_type column = entry->column();
1565
1566 // end of row reached? this should not happen, since we only
1567 // operate on compressed matrices!
1568 Assert(column != SparsityPattern::invalid_entry,
1569 ExcMatrixNotClosed());
1570
1571 if (distribute[column] == numbers::invalid_size_type)
1572 // distribute entry at irregular row @p row and regular
1573 // column column. set old entry to zero
1574 {
1575 for (size_type q = 0;
1576 q != lines[distribute[row]].entries.size();
1577 ++q)
1578 {
1579 // need a temporary variable to avoid errors like
1580 // no known conversion from 'complex<typename
1581 // ProductType<float, double>::type>' to 'const
1582 // complex<float>' for 3rd argument
1583 number v = static_cast<number>(entry->value());
1584 v *= lines[distribute[row]].entries[q].second;
1585 uncondensed.add(lines[distribute[row]].entries[q].first,
1586 column,
1587 v);
1588 }
1589
1590 // set old entry to zero
1591 entry->value() = 0.;
1592 }
1593 else
1594 // distribute entry at irregular row @p row and irregular
1595 // column @p column set old entry to one on main diagonal,
1596 // zero otherwise
1597 {
1598 for (size_type p = 0;
1599 p != lines[distribute[row]].entries.size();
1600 ++p)
1601 {
1602 for (size_type q = 0;
1603 q != lines[distribute[column]].entries.size();
1604 ++q)
1605 {
1606 // need a temporary variable to avoid errors like
1607 // no known conversion from 'complex<typename
1608 // ProductType<float, double>::type>' to 'const
1609 // complex<float>' for 3rd argument
1610 number v = static_cast<number>(entry->value());
1611 v *= lines[distribute[row]].entries[p].second *
1612 lines[distribute[column]].entries[q].second;
1613 uncondensed.add(
1614 lines[distribute[row]].entries[p].first,
1615 lines[distribute[column]].entries[q].first,
1616 v);
1617 }
1618
1619 if (use_vectors == true)
1620 vec(lines[distribute[row]].entries[p].first) -=
1621 static_cast<number>(entry->value()) *
1622 lines[distribute[row]].entries[p].second *
1623 lines[distribute[column]].inhomogeneity;
1624 }
1625
1626 // set old entry to correct value
1627 entry->value() = (row == column ? average_diagonal : 0.);
1628 }
1629 }
1630
1631 // take care of vector
1632 if (use_vectors == true)
1633 {
1634 for (size_type q = 0; q != lines[distribute[row]].entries.size();
1635 ++q)
1636 vec(lines[distribute[row]].entries[q].first) +=
1637 (vec(row) * lines[distribute[row]].entries[q].second);
1638
1639 vec(lines[distribute[row]].index) = 0.;
1640 }
1641 }
1642 }
1643 }
1644
1645
1646
1647 template <typename number>
1648 template <class BlockVectorType>
1649 void
condense(BlockSparseMatrix<number> & uncondensed,BlockVectorType & vec)1650 AffineConstraints<number>::condense(BlockSparseMatrix<number> &uncondensed,
1651 BlockVectorType & vec) const
1652 {
1653 // check whether we work on real vectors or we just used a dummy when
1654 // calling the other function above.
1655 const bool use_vectors = vec.n_blocks() == 0 ? false : true;
1656
1657 const size_type blocks = uncondensed.n_block_rows();
1658
1659 const BlockSparsityPattern &sparsity = uncondensed.get_sparsity_pattern();
1660
1661 Assert(sorted == true, ExcMatrixNotClosed());
1662 Assert(sparsity.is_compressed() == true, ExcMatrixNotClosed());
1663 Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
1664 Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
1665 Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
1666 Assert(sparsity.get_column_indices() == sparsity.get_row_indices(),
1667 ExcNotQuadratic());
1668
1669 if (use_vectors == true)
1670 {
1671 AssertDimension(vec.size(), sparsity.n_rows());
1672 AssertDimension(vec.n_blocks(), sparsity.n_block_rows());
1673 }
1674
1675 number average_diagonal = 0.;
1676 for (size_type b = 0; b < uncondensed.n_block_rows(); ++b)
1677 for (size_type i = 0; i < uncondensed.block(b, b).m(); ++i)
1678 average_diagonal += std::abs(uncondensed.block(b, b).diag_element(i));
1679 average_diagonal /= uncondensed.m();
1680
1681 const BlockIndices &index_mapping = sparsity.get_column_indices();
1682
1683 // store for each index whether it must be distributed or not. If entry
1684 // is numbers::invalid_size_type, no distribution is necessary.
1685 // otherwise, the number states which line in the constraints object
1686 // handles this index
1687 std::vector<size_type> distribute(sparsity.n_rows(),
1688 numbers::invalid_size_type);
1689
1690 for (size_type c = 0; c < lines.size(); ++c)
1691 distribute[lines[c].index] = c;
1692
1693 const size_type n_rows = sparsity.n_rows();
1694 for (size_type row = 0; row < n_rows; ++row)
1695 {
1696 // get index of this row within the blocks
1697 const std::pair<size_type, size_type> block_index =
1698 index_mapping.global_to_local(row);
1699 const size_type block_row = block_index.first;
1700
1701 if (distribute[row] == numbers::invalid_size_type)
1702 // regular line. loop over all columns and see whether this column
1703 // must be distributed
1704 {
1705 // to loop over all entries in this row, we have to loop over all
1706 // blocks in this blockrow and the corresponding row therein
1707 for (size_type block_col = 0; block_col < blocks; ++block_col)
1708 {
1709 for (typename SparseMatrix<number>::iterator entry =
1710 uncondensed.block(block_row, block_col)
1711 .begin(block_index.second);
1712 entry != uncondensed.block(block_row, block_col)
1713 .end(block_index.second);
1714 ++entry)
1715 {
1716 const size_type global_col =
1717 index_mapping.local_to_global(block_col, entry->column());
1718
1719 if (distribute[global_col] != numbers::invalid_size_type)
1720 // distribute entry at regular row @p row and irregular
1721 // column global_col; set old entry to zero
1722 {
1723 const number old_value = entry->value();
1724
1725 for (size_type q = 0;
1726 q != lines[distribute[global_col]].entries.size();
1727 ++q)
1728 uncondensed.add(
1729 row,
1730 lines[distribute[global_col]].entries[q].first,
1731 old_value *
1732 lines[distribute[global_col]].entries[q].second);
1733
1734 // need to subtract this element from the vector.
1735 // this corresponds to an explicit elimination in the
1736 // respective row of the inhomogeneous constraint in
1737 // the matrix with Gauss elimination
1738 if (use_vectors == true)
1739 vec(row) -= static_cast<number>(entry->value()) *
1740 lines[distribute[global_col]].inhomogeneity;
1741
1742 entry->value() = 0.;
1743 }
1744 }
1745 }
1746 }
1747 else
1748 {
1749 // row must be distributed. split the whole row into the chunks
1750 // defined by the blocks
1751 for (size_type block_col = 0; block_col < blocks; ++block_col)
1752 {
1753 for (typename SparseMatrix<number>::iterator entry =
1754 uncondensed.block(block_row, block_col)
1755 .begin(block_index.second);
1756 entry != uncondensed.block(block_row, block_col)
1757 .end(block_index.second);
1758 ++entry)
1759 {
1760 const size_type global_col =
1761 index_mapping.local_to_global(block_col, entry->column());
1762
1763 if (distribute[global_col] == numbers::invalid_size_type)
1764 // distribute entry at irregular row @p row and regular
1765 // column global_col. set old entry to zero
1766 {
1767 const number old_value = entry->value();
1768
1769 for (size_type q = 0;
1770 q != lines[distribute[row]].entries.size();
1771 ++q)
1772 uncondensed.add(
1773 lines[distribute[row]].entries[q].first,
1774 global_col,
1775 old_value * lines[distribute[row]].entries[q].second);
1776
1777 entry->value() = 0.;
1778 }
1779 else
1780 // distribute entry at irregular row @p row and
1781 // irregular column @p global_col set old entry to one
1782 // if on main diagonal, zero otherwise
1783 {
1784 const number old_value = entry->value();
1785
1786 for (size_type p = 0;
1787 p != lines[distribute[row]].entries.size();
1788 ++p)
1789 {
1790 for (size_type q = 0;
1791 q !=
1792 lines[distribute[global_col]].entries.size();
1793 ++q)
1794 uncondensed.add(
1795 lines[distribute[row]].entries[p].first,
1796 lines[distribute[global_col]].entries[q].first,
1797 old_value *
1798 lines[distribute[row]].entries[p].second *
1799 lines[distribute[global_col]]
1800 .entries[q]
1801 .second);
1802
1803 if (use_vectors == true)
1804 vec(lines[distribute[row]].entries[p].first) -=
1805 old_value *
1806 lines[distribute[row]].entries[p].second *
1807 lines[distribute[global_col]].inhomogeneity;
1808 }
1809
1810 entry->value() =
1811 (row == global_col ? average_diagonal : 0.);
1812 }
1813 }
1814 }
1815
1816 // take care of vector
1817 if (use_vectors == true)
1818 {
1819 for (size_type q = 0; q != lines[distribute[row]].entries.size();
1820 ++q)
1821 vec(lines[distribute[row]].entries[q].first) +=
1822 (vec(row) * lines[distribute[row]].entries[q].second);
1823
1824 vec(lines[distribute[row]].index) = 0.;
1825 }
1826 }
1827 }
1828 }
1829
1830
1831
1832 // TODO: I'm sure the following could be made more elegant by using a bit of
1833 // introspection using static member variables of the various vector
1834 // classes to dispatch between the different functions, rather than using
1835 // knowledge of the individual types
1836
1837 // number of functions to select the right implementation for set_zero().
1838 namespace internal
1839 {
1840 namespace AffineConstraintsImplementation
1841 {
1842 using size_type = types::global_dof_index;
1843
1844 template <class VectorType>
1845 void
1846 set_zero_parallel(const std::vector<size_type> &cm,
1847 VectorType & vec,
1848 size_type shift = 0)
1849 {
1850 Assert(!vec.has_ghost_elements(), ExcInternalError());
1851 IndexSet locally_owned = vec.locally_owned_elements();
1852 for (const size_type index : cm)
1853 {
1854 // If shift>0 then we are working on a part of a BlockVector
1855 // so vec(i) is actually the global entry i+shift.
1856 // We first make sure the line falls into the range of vec,
1857 // then check if is part of the local part of the vector, before
1858 // finally setting value to 0.
1859 if (index < shift)
1860 continue;
1861 const size_type idx = index - shift;
1862 if (idx < vec.size() && locally_owned.is_element(idx))
1863 internal::ElementAccess<VectorType>::set(0., idx, vec);
1864 }
1865 }
1866
1867 template <typename number>
1868 void
1869 set_zero_parallel(const std::vector<size_type> & cm,
1870 LinearAlgebra::distributed::Vector<number> &vec,
1871 size_type shift = 0)
1872 {
1873 for (const size_type index : cm)
1874 {
1875 // If shift>0 then we are working on a part of a BlockVector
1876 // so vec(i) is actually the global entry i+shift.
1877 // We first make sure the line falls into the range of vec,
1878 // then check if is part of the local part of the vector, before
1879 // finally setting the value to 0.
1880 if (index < shift)
1881 continue;
1882 const size_type idx = index - shift;
1883 if (vec.in_local_range(idx))
1884 vec(idx) = 0.;
1885 }
1886 vec.zero_out_ghosts();
1887 }
1888
1889 #ifdef DEAL_II_COMPILER_CUDA_AWARE
1890 template <typename Number>
1891 __global__ void
set_zero_kernel(const size_type * constrained_dofs,const unsigned int n_constrained_dofs,Number * dst)1892 set_zero_kernel(const size_type * constrained_dofs,
1893 const unsigned int n_constrained_dofs,
1894 Number * dst)
1895 {
1896 const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
1897 if (index < n_constrained_dofs)
1898 dst[constrained_dofs[index]] = 0;
1899 }
1900
1901 template <typename number>
1902 void
1903 set_zero_parallel(
1904 const std::vector<size_type> & cm,
1905 LinearAlgebra::distributed::Vector<number, MemorySpace::CUDA> &vec,
1906 size_type shift = 0)
1907 {
1908 Assert(shift == 0, ExcNotImplemented());
1909 (void)shift;
1910 std::vector<size_type> constrained_local_dofs_host;
1911 constrained_local_dofs_host.reserve(cm.size());
1912
1913 for (const auto global_index : cm)
1914 if (vec.in_local_range(global_index))
1915 constrained_local_dofs_host.push_back(
1916 vec.get_partitioner()->global_to_local(global_index));
1917
1918 const int n_constraints = constrained_local_dofs_host.size();
1919 size_type *constrained_local_dofs_device;
1920 Utilities::CUDA::malloc(constrained_local_dofs_device, n_constraints);
1921 Utilities::CUDA::copy_to_dev(constrained_local_dofs_host,
1922 constrained_local_dofs_device);
1923
1924 const int n_blocks = 1 + (n_constraints - 1) / CUDAWrappers::block_size;
1925 set_zero_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
1926 constrained_local_dofs_device, n_constraints, vec.get_values());
1927 AssertCudaKernel();
1928
1929 Utilities::CUDA::free(constrained_local_dofs_device);
1930
1931 vec.zero_out_ghosts();
1932 }
1933 #endif
1934
1935 template <class VectorType>
1936 void
set_zero_in_parallel(const std::vector<size_type> & cm,VectorType & vec,std::integral_constant<bool,false>)1937 set_zero_in_parallel(const std::vector<size_type> &cm,
1938 VectorType & vec,
1939 std::integral_constant<bool, false>)
1940 {
1941 set_zero_parallel(cm, vec, 0);
1942 }
1943
1944 // in parallel for BlockVectors
1945 template <class VectorType>
1946 void
set_zero_in_parallel(const std::vector<size_type> & cm,VectorType & vec,std::integral_constant<bool,true>)1947 set_zero_in_parallel(const std::vector<size_type> &cm,
1948 VectorType & vec,
1949 std::integral_constant<bool, true>)
1950 {
1951 size_type start_shift = 0;
1952 for (size_type j = 0; j < vec.n_blocks(); ++j)
1953 {
1954 set_zero_parallel(cm, vec.block(j), start_shift);
1955 start_shift += vec.block(j).size();
1956 }
1957 }
1958
1959 template <class VectorType>
1960 void
set_zero_serial(const std::vector<size_type> & cm,VectorType & vec)1961 set_zero_serial(const std::vector<size_type> &cm, VectorType &vec)
1962 {
1963 for (const size_type index : cm)
1964 vec(index) = 0.;
1965 }
1966
1967 template <class VectorType>
1968 void
set_zero_all(const std::vector<size_type> & cm,VectorType & vec)1969 set_zero_all(const std::vector<size_type> &cm, VectorType &vec)
1970 {
1971 set_zero_in_parallel<VectorType>(
1972 cm,
1973 vec,
1974 std::integral_constant<bool, IsBlockVector<VectorType>::value>());
1975 vec.compress(VectorOperation::insert);
1976 }
1977
1978 template <class T>
1979 void
set_zero_all(const std::vector<size_type> & cm,dealii::Vector<T> & vec)1980 set_zero_all(const std::vector<size_type> &cm, dealii::Vector<T> &vec)
1981 {
1982 set_zero_serial(cm, vec);
1983 }
1984
1985 template <class T>
1986 void
set_zero_all(const std::vector<size_type> & cm,dealii::BlockVector<T> & vec)1987 set_zero_all(const std::vector<size_type> &cm, dealii::BlockVector<T> &vec)
1988 {
1989 set_zero_serial(cm, vec);
1990 }
1991 } // namespace AffineConstraintsImplementation
1992 } // namespace internal
1993
1994 template <typename number>
1995 template <typename VectorType>
1996 void
distribute_local_to_global(const Vector<number> & local_vector,const std::vector<size_type> & local_dof_indices,VectorType & global_vector,const FullMatrix<number> & local_matrix)1997 AffineConstraints<number>::distribute_local_to_global(
1998 const Vector<number> & local_vector,
1999 const std::vector<size_type> &local_dof_indices,
2000 VectorType & global_vector,
2001 const FullMatrix<number> & local_matrix) const
2002 {
2003 distribute_local_to_global(local_vector,
2004 local_dof_indices,
2005 local_dof_indices,
2006 global_vector,
2007 local_matrix,
2008 true);
2009 }
2010
2011 template <typename number>
2012 template <typename VectorType>
2013 void
distribute_local_to_global(const Vector<number> & local_vector,const std::vector<size_type> & local_dof_indices_row,const std::vector<size_type> & local_dof_indices_col,VectorType & global_vector,const FullMatrix<number> & local_matrix,bool diagonal)2014 AffineConstraints<number>::distribute_local_to_global(
2015 const Vector<number> & local_vector,
2016 const std::vector<size_type> &local_dof_indices_row,
2017 const std::vector<size_type> &local_dof_indices_col,
2018 VectorType & global_vector,
2019 const FullMatrix<number> & local_matrix,
2020 bool diagonal) const
2021 {
2022 Assert(sorted == true, ExcMatrixNotClosed());
2023 AssertDimension(local_vector.size(), local_dof_indices_row.size());
2024 AssertDimension(local_matrix.m(), local_dof_indices_row.size());
2025 AssertDimension(local_matrix.n(), local_dof_indices_col.size());
2026
2027 // diagonal checks if we have only one index set (if both are equal
2028 // diagonal should be set to true).
2029 // If true we do both, assembly of the right hand side (next lines)
2030 // and (see further below) modifications of the right hand side
2031 // according to the inhomogeneous constraints.
2032 // Otherwise we only modify the right hand side according to
2033 // local_matrix and the inhomogeneous constraints, and omit the vector add.
2034
2035 const size_type m_local_dofs = local_dof_indices_row.size();
2036 const size_type n_local_dofs = local_dof_indices_col.size();
2037 if (lines.empty())
2038 {
2039 if (diagonal)
2040 global_vector.add(local_dof_indices_row, local_vector);
2041 }
2042 else
2043 for (size_type i = 0; i < n_local_dofs; ++i)
2044 {
2045 // check whether the current index is
2046 // constrained. if not, just write the entry
2047 // into the vector. otherwise, need to resolve
2048 // the constraint
2049 if (is_constrained(local_dof_indices_col[i]) == false)
2050 {
2051 if (diagonal)
2052 global_vector(local_dof_indices_row[i]) += local_vector(i);
2053 continue;
2054 }
2055
2056 // find the constraint line to the given
2057 // global dof index
2058 const size_type line_index =
2059 calculate_line_index(local_dof_indices_col[i]);
2060 AssertIndexRange(line_index, lines_cache.size());
2061 AssertIndexRange(lines_cache[line_index], lines.size());
2062 const ConstraintLine &position = lines[lines_cache[line_index]];
2063
2064 // Gauss elimination of the matrix columns with the inhomogeneity.
2065 // Go through them one by one and again check whether they are
2066 // constrained. If so, distribute the constraint
2067 const auto val = position.inhomogeneity;
2068 if (val != number(0.))
2069 for (size_type j = 0; j < m_local_dofs; ++j)
2070 {
2071 if (is_constrained(local_dof_indices_row[j]) == false)
2072 {
2073 global_vector(local_dof_indices_row[j]) -=
2074 val * local_matrix(j, i);
2075 continue;
2076 }
2077
2078 const number matrix_entry = local_matrix(j, i);
2079
2080 if (matrix_entry == number())
2081 continue;
2082
2083 const ConstraintLine &position_j =
2084 lines[lines_cache[calculate_line_index(
2085 local_dof_indices_row[j])]];
2086
2087 for (size_type q = 0; q < position_j.entries.size(); ++q)
2088 {
2089 Assert(!(!local_lines.size() ||
2090 local_lines.is_element(
2091 position_j.entries[q].first)) ||
2092 is_constrained(position_j.entries[q].first) == false,
2093 ExcMessage("Tried to distribute to a fixed dof."));
2094 global_vector(position_j.entries[q].first) -=
2095 val * position_j.entries[q].second * matrix_entry;
2096 }
2097 }
2098
2099 // now distribute the constraint,
2100 // but make sure we don't touch
2101 // the entries of fixed dofs
2102 if (diagonal)
2103 {
2104 for (size_type j = 0; j < position.entries.size(); ++j)
2105 {
2106 Assert(!(!local_lines.size() ||
2107 local_lines.is_element(position.entries[j].first)) ||
2108 is_constrained(position.entries[j].first) == false,
2109 ExcMessage("Tried to distribute to a fixed dof."));
2110 global_vector(position.entries[j].first) +=
2111 local_vector(i) * position.entries[j].second;
2112 }
2113 }
2114 }
2115 }
2116
2117 namespace internal
2118 {
2119 // create an output vector that consists of the input vector's locally owned
2120 // elements plus some ghost elements that need to be imported from elsewhere
2121 //
2122 // this is an operation that is different for all vector types and so we
2123 // need a few overloads
2124 #ifdef DEAL_II_WITH_TRILINOS
2125 inline void
import_vector_with_ghost_elements(const TrilinosWrappers::MPI::Vector & vec,const IndexSet &,const IndexSet & needed_elements,TrilinosWrappers::MPI::Vector & output,const std::integral_constant<bool,false>)2126 import_vector_with_ghost_elements(
2127 const TrilinosWrappers::MPI::Vector &vec,
2128 const IndexSet & /*locally_owned_elements*/,
2129 const IndexSet & needed_elements,
2130 TrilinosWrappers::MPI::Vector &output,
2131 const std::integral_constant<bool, false> /*is_block_vector*/)
2132 {
2133 Assert(!vec.has_ghost_elements(), ExcGhostsPresent());
2134 # ifdef DEAL_II_WITH_MPI
2135 const Epetra_MpiComm *mpi_comm =
2136 dynamic_cast<const Epetra_MpiComm *>(&vec.trilinos_vector().Comm());
2137
2138 Assert(mpi_comm != nullptr, ExcInternalError());
2139 output.reinit(needed_elements, mpi_comm->GetMpiComm());
2140 # else
2141 output.reinit(needed_elements, MPI_COMM_SELF);
2142 # endif
2143 output = vec;
2144 }
2145 #endif
2146
2147 #ifdef DEAL_II_WITH_PETSC
2148 inline void
import_vector_with_ghost_elements(const PETScWrappers::MPI::Vector & vec,const IndexSet & locally_owned_elements,const IndexSet & needed_elements,PETScWrappers::MPI::Vector & output,const std::integral_constant<bool,false>)2149 import_vector_with_ghost_elements(
2150 const PETScWrappers::MPI::Vector &vec,
2151 const IndexSet & locally_owned_elements,
2152 const IndexSet & needed_elements,
2153 PETScWrappers::MPI::Vector & output,
2154 const std::integral_constant<bool, false> /*is_block_vector*/)
2155 {
2156 output.reinit(locally_owned_elements,
2157 needed_elements,
2158 vec.get_mpi_communicator());
2159 output = vec;
2160 }
2161 #endif
2162
2163 template <typename number>
2164 void
import_vector_with_ghost_elements(const LinearAlgebra::distributed::Vector<number> & vec,const IndexSet & locally_owned_elements,const IndexSet & needed_elements,LinearAlgebra::distributed::Vector<number> & output,const std::integral_constant<bool,false>)2165 import_vector_with_ghost_elements(
2166 const LinearAlgebra::distributed::Vector<number> &vec,
2167 const IndexSet & locally_owned_elements,
2168 const IndexSet & needed_elements,
2169 LinearAlgebra::distributed::Vector<number> & output,
2170 const std::integral_constant<bool, false> /*is_block_vector*/)
2171 {
2172 // TODO: the in vector might already have all elements. need to find a
2173 // way to efficiently avoid the copy then
2174 const_cast<LinearAlgebra::distributed::Vector<number> &>(vec)
2175 .zero_out_ghosts();
2176 output.reinit(locally_owned_elements,
2177 needed_elements,
2178 vec.get_mpi_communicator());
2179 output = vec;
2180 output.update_ghost_values();
2181 }
2182
2183 // all other vector non-block vector types are sequential and we should
2184 // not have this function called at all -- so throw an exception
2185 template <typename Vector>
2186 void
import_vector_with_ghost_elements(const Vector &,const IndexSet &,const IndexSet &,Vector &,const std::integral_constant<bool,false>)2187 import_vector_with_ghost_elements(
2188 const Vector & /*vec*/,
2189 const IndexSet & /*locally_owned_elements*/,
2190 const IndexSet & /*needed_elements*/,
2191 Vector & /*output*/,
2192 const std::integral_constant<bool, false> /*is_block_vector*/)
2193 {
2194 Assert(false, ExcMessage("We shouldn't even get here!"));
2195 }
2196
2197 // for block vectors, simply dispatch to the individual blocks
2198 template <class VectorType>
2199 void
import_vector_with_ghost_elements(const VectorType & vec,const IndexSet & locally_owned_elements,const IndexSet & needed_elements,VectorType & output,const std::integral_constant<bool,true>)2200 import_vector_with_ghost_elements(
2201 const VectorType &vec,
2202 const IndexSet & locally_owned_elements,
2203 const IndexSet & needed_elements,
2204 VectorType & output,
2205 const std::integral_constant<bool, true> /*is_block_vector*/)
2206 {
2207 output.reinit(vec.n_blocks());
2208
2209 types::global_dof_index block_start = 0;
2210 for (unsigned int b = 0; b < vec.n_blocks(); ++b)
2211 {
2212 import_vector_with_ghost_elements(
2213 vec.block(b),
2214 locally_owned_elements.get_view(block_start,
2215 block_start + vec.block(b).size()),
2216 needed_elements.get_view(block_start,
2217 block_start + vec.block(b).size()),
2218 output.block(b),
2219 std::integral_constant<bool, false>());
2220 block_start += vec.block(b).size();
2221 }
2222
2223 output.collect_sizes();
2224 }
2225 } // namespace internal
2226
2227 template <typename number>
2228 template <class VectorType>
2229 void
distribute(VectorType & vec)2230 AffineConstraints<number>::distribute(VectorType &vec) const
2231 {
2232 Assert(sorted == true, ExcMatrixNotClosed());
2233
2234 // if the vector type supports parallel storage and if the vector actually
2235 // does store only part of the vector, distributing is slightly more
2236 // complicated. we might be able to skip the complicated part if one
2237 // processor owns everything and pretend that this is a sequential vector,
2238 // but it is difficult for the other processors to know whether they should
2239 // not do anything or if other processors will create a temporary vector,
2240 // exchange data (requiring communication, maybe even with the processors
2241 // that do not own anything because of that particular parallel model), and
2242 // call compress() finally. the first case here is for the complicated case,
2243 // the last else is for the simple case (sequential vector)
2244 const IndexSet vec_owned_elements = vec.locally_owned_elements();
2245
2246 if (dealii::is_serial_vector<VectorType>::value == false)
2247 {
2248 // This processor owns only part of the vector. one may think that
2249 // every processor should be able to simply communicate those elements
2250 // it owns and for which it knows that they act as sources to constrained
2251 // DoFs to the owner of these DoFs. This would lead to a scheme where all
2252 // we need to do is to add some local elements to (possibly non-local)
2253 // ones and then call compress().
2254 //
2255 // Alas, this scheme does not work as evidenced by the disaster of bug
2256 // #51, see http://code.google.com/p/dealii/issues/detail?id=51 and the
2257 // reversion of one attempt that implements this in r29662. Rather, we
2258 // need to get a vector that has all the *sources* or constraints we
2259 // own locally, possibly as ghost vector elements, then read from them,
2260 // and finally throw away the ghosted vector. Implement this in the
2261 // following.
2262 IndexSet needed_elements = vec_owned_elements;
2263
2264 for (const ConstraintLine &line : lines)
2265 if (vec_owned_elements.is_element(line.index))
2266 for (const std::pair<size_type, number> &entry : line.entries)
2267 if (!vec_owned_elements.is_element(entry.first))
2268 needed_elements.add_index(entry.first);
2269
2270 VectorType ghosted_vector;
2271 internal::import_vector_with_ghost_elements(
2272 vec,
2273 vec_owned_elements,
2274 needed_elements,
2275 ghosted_vector,
2276 std::integral_constant<bool, IsBlockVector<VectorType>::value>());
2277
2278 for (const ConstraintLine &line : lines)
2279 if (vec_owned_elements.is_element(line.index))
2280 {
2281 typename VectorType::value_type new_value = line.inhomogeneity;
2282 for (const std::pair<size_type, number> &entry : line.entries)
2283 new_value +=
2284 (static_cast<typename VectorType::value_type>(
2285 internal::ElementAccess<VectorType>::get(ghosted_vector,
2286 entry.first)) *
2287 entry.second);
2288 AssertIsFinite(new_value);
2289 internal::ElementAccess<VectorType>::set(new_value,
2290 line.index,
2291 vec);
2292 }
2293
2294 // now compress to communicate the entries that we added to
2295 // and that weren't to local processors to the owner
2296 //
2297 // this shouldn't be strictly necessary but it probably doesn't
2298 // hurt either
2299 vec.compress(VectorOperation::insert);
2300 }
2301 else
2302 // purely sequential vector (either because the type doesn't
2303 // support anything else or because it's completely stored
2304 // locally)
2305 {
2306 for (const ConstraintLine &next_constraint : lines)
2307 {
2308 // fill entry in line
2309 // next_constraint.index by adding the
2310 // different contributions
2311 typename VectorType::value_type new_value =
2312 next_constraint.inhomogeneity;
2313 for (const std::pair<size_type, number> &entry :
2314 next_constraint.entries)
2315 new_value +=
2316 (static_cast<typename VectorType::value_type>(
2317 internal::ElementAccess<VectorType>::get(vec, entry.first)) *
2318 entry.second);
2319 AssertIsFinite(new_value);
2320 internal::ElementAccess<VectorType>::set(new_value,
2321 next_constraint.index,
2322 vec);
2323 }
2324 }
2325 }
2326
2327 // Some helper definitions for the local_to_global functions.
2328 namespace internal
2329 {
2330 namespace AffineConstraints
2331 {
Distributing(const size_type global_row,const size_type local_row)2332 inline Distributing::Distributing(const size_type global_row,
2333 const size_type local_row)
2334 : global_row(global_row)
2335 , local_row(local_row)
2336 , constraint_position(numbers::invalid_size_type)
2337 {}
2338
2339
2340
Distributing(const Distributing & in)2341 inline Distributing::Distributing(const Distributing &in)
2342 : constraint_position(numbers::invalid_size_type)
2343 {
2344 *this = (in);
2345 }
2346
2347
2348
2349 inline Distributing &
2350 Distributing::operator=(const Distributing &in)
2351 {
2352 global_row = in.global_row;
2353 local_row = in.local_row;
2354 // the constraints pointer should not contain any data here.
2355 Assert(constraint_position == numbers::invalid_size_type,
2356 ExcInternalError());
2357
2358 if (in.constraint_position != numbers::invalid_size_type)
2359 {
2360 constraint_position = in.constraint_position;
2361 in.constraint_position = numbers::invalid_size_type;
2362 }
2363 return *this;
2364 }
2365
2366
2367
2368 template <typename number>
DataCache()2369 DataCache<number>::DataCache()
2370 : row_length(8)
2371 {}
2372
2373
2374
2375 template <typename number>
2376 void
reinit()2377 DataCache<number>::reinit()
2378 {
2379 individual_size.resize(0);
2380 data.resize(0);
2381 }
2382
2383
2384
2385 template <typename number>
2386 size_type
insert_new_index(const std::pair<size_type,number> & pair)2387 DataCache<number>::insert_new_index(
2388 const std::pair<size_type, number> &pair)
2389 {
2390 Assert(row_length > 0, ExcInternalError());
2391 const unsigned int index = individual_size.size();
2392 individual_size.push_back(1);
2393 data.resize(individual_size.size() * row_length);
2394 data[index * row_length] = pair;
2395 individual_size[index] = 1;
2396 return index;
2397 }
2398
2399
2400
2401 template <typename number>
2402 void
append_index(const size_type index,const std::pair<size_type,number> & pair)2403 DataCache<number>::append_index(const size_type index,
2404 const std::pair<size_type, number> &pair)
2405 {
2406 AssertIndexRange(index, individual_size.size());
2407 const size_type my_length = individual_size[index];
2408 if (my_length == row_length)
2409 {
2410 AssertDimension(data.size(), individual_size.size() * row_length);
2411 // no space left in this row, need to double row_length and
2412 // rearrange the data items. Move all items to the right except the
2413 // first one, starting at the back. Since individual_size contains
2414 // at least one element when we get here, subtracting 1 works fine.
2415 data.resize(2 * data.size());
2416 for (size_type i = individual_size.size() - 1; i > 0; --i)
2417 {
2418 const auto ptr = data.data();
2419 std::move_backward(ptr + i * row_length,
2420 ptr + i * row_length + individual_size[i],
2421 ptr + i * 2 * row_length + individual_size[i]);
2422 }
2423 row_length *= 2;
2424 }
2425 data[index * row_length + my_length] = pair;
2426 individual_size[index] = my_length + 1;
2427 }
2428
2429
2430
2431 template <typename number>
2432 size_type
get_size(const size_type index)2433 DataCache<number>::get_size(const size_type index) const
2434 {
2435 return individual_size[index];
2436 }
2437
2438
2439
2440 template <typename number>
2441 const std::pair<size_type, number> *
get_entry(const size_type index)2442 DataCache<number>::get_entry(const size_type index) const
2443 {
2444 return &data[index * row_length];
2445 }
2446
2447
2448
2449 template <typename number>
GlobalRowsFromLocal()2450 GlobalRowsFromLocal<number>::GlobalRowsFromLocal()
2451 : n_active_rows(0)
2452 , n_inhomogeneous_rows(0)
2453 {}
2454
2455
2456
2457 template <typename number>
2458 void
reinit(const size_type n_local_rows)2459 GlobalRowsFromLocal<number>::reinit(const size_type n_local_rows)
2460 {
2461 total_row_indices.resize(n_local_rows);
2462 for (unsigned int i = 0; i < n_local_rows; ++i)
2463 total_row_indices[i].constraint_position = numbers::invalid_size_type;
2464 n_active_rows = n_local_rows;
2465 n_inhomogeneous_rows = 0;
2466 data_cache.reinit();
2467 }
2468
2469
2470
2471 template <typename number>
2472 void
print(std::ostream & os)2473 GlobalRowsFromLocal<number>::print(std::ostream &os)
2474 {
2475 os << "Active rows " << n_active_rows << std::endl
2476 << "Constr rows " << n_constraints() << std::endl
2477 << "Inhom rows " << n_inhomogeneous_rows << std::endl
2478 << "Local: ";
2479 for (const auto &total_row_index : total_row_indices)
2480 os << ' ' << std::setw(4) << total_row_index.local_row;
2481 os << std::endl << "Global:";
2482 for (const auto &total_row_index : total_row_indices)
2483 os << ' ' << std::setw(4) << total_row_index.global_row;
2484 os << std::endl << "ConPos:";
2485 for (const auto &total_row_index : total_row_indices)
2486 os << ' ' << std::setw(4) << total_row_index.constraint_position;
2487 os << std::endl;
2488 }
2489
2490
2491
2492 template <typename number>
2493 size_type
size()2494 GlobalRowsFromLocal<number>::size() const
2495 {
2496 return n_active_rows;
2497 }
2498
2499
2500
2501 template <typename number>
2502 size_type
size(const size_type counter_index)2503 GlobalRowsFromLocal<number>::size(const size_type counter_index) const
2504 {
2505 return (total_row_indices[counter_index].constraint_position ==
2506 numbers::invalid_size_type ?
2507 0 :
2508 data_cache.get_size(
2509 total_row_indices[counter_index].constraint_position));
2510 }
2511
2512
2513
2514 template <typename number>
2515 size_type
global_row(const size_type counter_index)2516 GlobalRowsFromLocal<number>::global_row(const size_type counter_index) const
2517 {
2518 return total_row_indices[counter_index].global_row;
2519 }
2520
2521
2522
2523 template <typename number>
2524 size_type &
global_row(const size_type counter_index)2525 GlobalRowsFromLocal<number>::global_row(const size_type counter_index)
2526 {
2527 return total_row_indices[counter_index].global_row;
2528 }
2529
2530
2531
2532 template <typename number>
2533 size_type
local_row(const size_type counter_index)2534 GlobalRowsFromLocal<number>::local_row(const size_type counter_index) const
2535 {
2536 return total_row_indices[counter_index].local_row;
2537 }
2538
2539
2540
2541 template <typename number>
2542 size_type &
local_row(const size_type counter_index)2543 GlobalRowsFromLocal<number>::local_row(const size_type counter_index)
2544 {
2545 return total_row_indices[counter_index].local_row;
2546 }
2547
2548
2549
2550 template <typename number>
2551 size_type
local_row(const size_type counter_index,const size_type index_in_constraint)2552 GlobalRowsFromLocal<number>::local_row(
2553 const size_type counter_index,
2554 const size_type index_in_constraint) const
2555 {
2556 return (data_cache.get_entry(total_row_indices[counter_index]
2557 .constraint_position)[index_in_constraint])
2558 .first;
2559 }
2560
2561
2562
2563 template <typename number>
2564 number
constraint_value(const size_type counter_index,const size_type index_in_constraint)2565 GlobalRowsFromLocal<number>::constraint_value(
2566 const size_type counter_index,
2567 const size_type index_in_constraint) const
2568 {
2569 return (data_cache.get_entry(total_row_indices[counter_index]
2570 .constraint_position)[index_in_constraint])
2571 .second;
2572 }
2573
2574
2575
2576 template <typename number>
2577 bool
have_indirect_rows()2578 GlobalRowsFromLocal<number>::have_indirect_rows() const
2579 {
2580 return data_cache.individual_size.empty() == false;
2581 }
2582
2583
2584
2585 template <typename number>
2586 void
insert_constraint(const size_type constrained_local_dof)2587 GlobalRowsFromLocal<number>::insert_constraint(
2588 const size_type constrained_local_dof)
2589 {
2590 --n_active_rows;
2591 total_row_indices[n_active_rows].local_row = constrained_local_dof;
2592 total_row_indices[n_active_rows].global_row = numbers::invalid_size_type;
2593 }
2594
2595
2596
2597 template <typename number>
2598 size_type
n_constraints()2599 GlobalRowsFromLocal<number>::n_constraints() const
2600 {
2601 return total_row_indices.size() - n_active_rows;
2602 }
2603
2604
2605
2606 template <typename number>
2607 size_type
n_inhomogeneities()2608 GlobalRowsFromLocal<number>::n_inhomogeneities() const
2609 {
2610 return n_inhomogeneous_rows;
2611 }
2612
2613
2614
2615 template <typename number>
2616 void
set_ith_constraint_inhomogeneous(const size_type i)2617 GlobalRowsFromLocal<number>::set_ith_constraint_inhomogeneous(
2618 const size_type i)
2619 {
2620 Assert(i >= n_inhomogeneous_rows, ExcInternalError());
2621 std::swap(total_row_indices[n_active_rows + i],
2622 total_row_indices[n_active_rows + n_inhomogeneous_rows]);
2623 n_inhomogeneous_rows++;
2624 }
2625
2626
2627
2628 template <typename number>
2629 size_type
constraint_origin(size_type i)2630 GlobalRowsFromLocal<number>::constraint_origin(size_type i) const
2631 {
2632 return total_row_indices[n_active_rows + i].local_row;
2633 }
2634
2635
2636 // a function that appends an additional row to the list of values, or
2637 // appends a value to an already existing row. Similar functionality as for
2638 // std::map<size_type,Distributing>, but here done for a
2639 // std::vector<Distributing>, much faster for short lists as we have them
2640 // here
2641 template <typename number>
2642 inline void
insert_index(const size_type global_row,const size_type local_row,const number constraint_value)2643 GlobalRowsFromLocal<number>::insert_index(const size_type global_row,
2644 const size_type local_row,
2645 const number constraint_value)
2646 {
2647 using index_iterator = std::vector<Distributing>::iterator;
2648 index_iterator pos, pos1;
2649 Distributing row_value(global_row);
2650 std::pair<size_type, number> constraint(local_row, constraint_value);
2651
2652 // check whether the list was really sorted before entering here
2653 for (size_type i = 1; i < n_active_rows; ++i)
2654 Assert(total_row_indices[i - 1] < total_row_indices[i],
2655 ExcInternalError());
2656
2657 pos = Utilities::lower_bound(total_row_indices.begin(),
2658 total_row_indices.begin() + n_active_rows,
2659 row_value);
2660 if (pos->global_row == global_row)
2661 pos1 = pos;
2662 else
2663 {
2664 pos1 = total_row_indices.insert(pos, row_value);
2665 ++n_active_rows;
2666 }
2667
2668 if (pos1->constraint_position == numbers::invalid_size_type)
2669 pos1->constraint_position = data_cache.insert_new_index(constraint);
2670 else
2671 data_cache.append_index(pos1->constraint_position, constraint);
2672 }
2673
2674 // this sort algorithm sorts std::vector<Distributing>, but does not take
2675 // the constraints into account. this means that in case that constraints
2676 // are already inserted, this function does not work as expected. Use
2677 // shellsort, which is very fast in case the indices are already sorted
2678 // (which is the usual case with DG elements), and not too slow in other
2679 // cases
2680 template <typename number>
2681 inline void
sort()2682 GlobalRowsFromLocal<number>::sort()
2683 {
2684 size_type i, j, j2, temp, templ, istep;
2685 size_type step;
2686
2687 // check whether the constraints are really empty.
2688 const size_type length = size();
2689
2690 // make sure that we are in the range of the vector
2691 AssertIndexRange(length, total_row_indices.size() + 1);
2692 for (size_type i = 0; i < length; ++i)
2693 Assert(total_row_indices[i].constraint_position ==
2694 numbers::invalid_size_type,
2695 ExcInternalError());
2696
2697 step = length / 2;
2698 while (step > 0)
2699 {
2700 for (i = step; i < length; i++)
2701 {
2702 istep = step;
2703 j = i;
2704 j2 = j - istep;
2705 temp = total_row_indices[i].global_row;
2706 templ = total_row_indices[i].local_row;
2707 if (total_row_indices[j2].global_row > temp)
2708 {
2709 while ((j >= istep) &&
2710 (total_row_indices[j2].global_row > temp))
2711 {
2712 total_row_indices[j].global_row =
2713 total_row_indices[j2].global_row;
2714 total_row_indices[j].local_row =
2715 total_row_indices[j2].local_row;
2716 j = j2;
2717 j2 -= istep;
2718 }
2719 total_row_indices[j].global_row = temp;
2720 total_row_indices[j].local_row = templ;
2721 }
2722 }
2723 step = step >> 1;
2724 }
2725 }
2726
2727
2728
2729 /**
2730 * This class is an accessor class to scratch data that is used
2731 * during calls to distribute_local_to_global and
2732 * add_entries_local_to_global. In order to avoid frequent memory
2733 * allocation, we keep the data alive from one call to the next in
2734 * a static variable. Since we want to allow for different number
2735 * types in matrices, this is a template.
2736 *
2737 * Since each thread gets its private version of scratch data out of the
2738 * ThreadLocalStorage, no conflicting access can occur. For this to be
2739 * valid, we need to make sure that no call within
2740 * distribute_local_to_global is made that by itself can spawn tasks.
2741 * Otherwise, we might end up in a situation where several threads fight for
2742 * the data.
2743 */
2744 template <typename number>
2745 class ScratchDataAccessor
2746 {
2747 public:
2748 /**
2749 * Constructor. Takes the scratch data object for the current
2750 * thread out of the provided object and marks it as used.
2751 */
ScratchDataAccessor(Threads::ThreadLocalStorage<ScratchData<number>> & tls_scratch_data)2752 ScratchDataAccessor(
2753 Threads::ThreadLocalStorage<ScratchData<number>> &tls_scratch_data)
2754 : my_scratch_data(&tls_scratch_data.get())
2755 {
2756 Assert(my_scratch_data->in_use == false,
2757 ExcMessage(
2758 "Access to thread-local scratch data tried, but it is already "
2759 "in use"));
2760 my_scratch_data->in_use = true;
2761 }
2762
2763 /**
2764 * Destructor. Mark scratch data as available again.
2765 */
~ScratchDataAccessor()2766 ~ScratchDataAccessor()
2767 {
2768 my_scratch_data->in_use = false;
2769 }
2770
2771 /**
2772 * Dereferencing operator.
2773 */
2774 ScratchData<number> &operator*()
2775 {
2776 return *my_scratch_data;
2777 }
2778
2779 /**
2780 * Dereferencing operator.
2781 */
2782 ScratchData<number> *operator->()
2783 {
2784 return my_scratch_data;
2785 }
2786
2787 private:
2788 ScratchData<number> *my_scratch_data;
2789 };
2790
2791
2792 // function for block matrices: Find out where in the list of local dofs
2793 // (sorted according to global ids) the individual blocks start.
2794 // Transform the global indices to block-local indices in order to be
2795 // able to use functions like vector.block(1)(block_local_id), instead of
2796 // vector(global_id). This avoids transforming indices one-by-one later
2797 // on.
2798 template <typename number, class BlockType>
2799 inline void
make_block_starts(const BlockType & block_object,GlobalRowsFromLocal<number> & global_rows,std::vector<size_type> & block_starts)2800 make_block_starts(const BlockType & block_object,
2801 GlobalRowsFromLocal<number> &global_rows,
2802 std::vector<size_type> & block_starts)
2803 {
2804 AssertDimension(block_starts.size(), block_object.n_block_rows() + 1);
2805
2806 using row_iterator = std::vector<Distributing>::iterator;
2807 row_iterator block_indices = global_rows.total_row_indices.begin();
2808
2809 const size_type num_blocks = block_object.n_block_rows();
2810 const size_type n_active_rows = global_rows.size();
2811
2812 // find end of rows.
2813 block_starts[0] = 0;
2814 for (size_type i = 1; i < num_blocks; ++i)
2815 {
2816 row_iterator first_block = Utilities::lower_bound(
2817 block_indices,
2818 global_rows.total_row_indices.begin() + n_active_rows,
2819 Distributing(block_object.get_row_indices().block_start(i)));
2820 block_starts[i] = first_block - global_rows.total_row_indices.begin();
2821 block_indices = first_block;
2822 }
2823 block_starts[num_blocks] = n_active_rows;
2824
2825 // transform row indices to block-local index space
2826 for (size_type i = block_starts[1]; i < n_active_rows; ++i)
2827 global_rows.global_row(i) =
2828 block_object.get_row_indices()
2829 .global_to_local(global_rows.global_row(i))
2830 .second;
2831 }
2832
2833 // same as before, but for std::vector<uint> instead of
2834 // GlobalRowsFromLocal. Used in functions for sparsity patterns.
2835 template <class BlockType>
2836 inline void
make_block_starts(const BlockType & block_object,std::vector<size_type> & row_indices,std::vector<size_type> & block_starts)2837 make_block_starts(const BlockType & block_object,
2838 std::vector<size_type> &row_indices,
2839 std::vector<size_type> &block_starts)
2840 {
2841 AssertDimension(block_starts.size(), block_object.n_block_rows() + 1);
2842
2843 using row_iterator = std::vector<size_type>::iterator;
2844 row_iterator col_indices = row_indices.begin();
2845
2846 const size_type num_blocks = block_object.n_block_rows();
2847
2848 // find end of rows.
2849 block_starts[0] = 0;
2850 for (size_type i = 1; i < num_blocks; ++i)
2851 {
2852 row_iterator first_block =
2853 Utilities::lower_bound(col_indices,
2854 row_indices.end(),
2855 block_object.get_row_indices().block_start(
2856 i));
2857 block_starts[i] = first_block - row_indices.begin();
2858 col_indices = first_block;
2859 }
2860 block_starts[num_blocks] = row_indices.size();
2861
2862 // transform row indices to local index space
2863 for (size_type i = block_starts[1]; i < row_indices.size(); ++i)
2864 row_indices[i] =
2865 block_object.get_row_indices().global_to_local(row_indices[i]).second;
2866 }
2867
2868 // resolves constraints of one column at the innermost loop. goes through
2869 // the origin of each global entry and finds out which data we need to
2870 // collect.
2871 template <typename number>
2872 static inline number
resolve_matrix_entry(const GlobalRowsFromLocal<number> & global_rows,const GlobalRowsFromLocal<number> & global_cols,const size_type i,const size_type j,const size_type loc_row,const FullMatrix<number> & local_matrix)2873 resolve_matrix_entry(const GlobalRowsFromLocal<number> &global_rows,
2874 const GlobalRowsFromLocal<number> &global_cols,
2875 const size_type i,
2876 const size_type j,
2877 const size_type loc_row,
2878 const FullMatrix<number> & local_matrix)
2879 {
2880 const size_type loc_col = global_cols.local_row(j);
2881 number col_val;
2882
2883 // case 1: row has direct contribution in local matrix. decide whether col
2884 // has a direct contribution. if not, set the value to zero.
2885 if (loc_row != numbers::invalid_size_type)
2886 {
2887 col_val = ((loc_col != numbers::invalid_size_type) ?
2888 local_matrix(loc_row, loc_col) :
2889 0);
2890
2891 // account for indirect contributions by constraints in column
2892 for (size_type p = 0; p < global_cols.size(j); ++p)
2893 col_val += (local_matrix(loc_row, global_cols.local_row(j, p)) *
2894 global_cols.constraint_value(j, p));
2895 }
2896
2897 // case 2: row has no direct contribution in local matrix
2898 else
2899 col_val = 0;
2900
2901 // account for indirect contributions by constraints in row, going trough
2902 // the direct and indirect references in the given column.
2903 for (size_type q = 0; q < global_rows.size(i); ++q)
2904 {
2905 number add_this =
2906 (loc_col != numbers::invalid_size_type) ?
2907 local_matrix(global_rows.local_row(i, q), loc_col) :
2908 0;
2909
2910 for (size_type p = 0; p < global_cols.size(j); ++p)
2911 add_this += (local_matrix(global_rows.local_row(i, q),
2912 global_cols.local_row(j, p)) *
2913 global_cols.constraint_value(j, p));
2914 col_val += add_this * global_rows.constraint_value(i, q);
2915 }
2916 return col_val;
2917 }
2918
2919 // computes all entries that need to be written into global_rows[i]. Lists
2920 // the resulting values in val_ptr, and the corresponding column indices in
2921 // col_ptr.
2922 template <typename number>
2923 inline void
resolve_matrix_row(const GlobalRowsFromLocal<number> & global_rows,const GlobalRowsFromLocal<number> & global_cols,const size_type i,const size_type column_start,const size_type column_end,const FullMatrix<number> & local_matrix,size_type * & col_ptr,number * & val_ptr)2924 resolve_matrix_row(const GlobalRowsFromLocal<number> &global_rows,
2925 const GlobalRowsFromLocal<number> &global_cols,
2926 const size_type i,
2927 const size_type column_start,
2928 const size_type column_end,
2929 const FullMatrix<number> & local_matrix,
2930 size_type *& col_ptr,
2931 number *& val_ptr)
2932 {
2933 if (column_end == column_start)
2934 return;
2935
2936 AssertIndexRange(column_end - 1, global_cols.size());
2937 const size_type loc_row = global_rows.local_row(i);
2938
2939 // fast function if there are no indirect references to any of the local
2940 // rows at all on this set of dofs (saves a lot of checks). the only check
2941 // we actually need to perform is whether the matrix element is zero.
2942 if (global_rows.have_indirect_rows() == false &&
2943 global_cols.have_indirect_rows() == false)
2944 {
2945 AssertIndexRange(loc_row, local_matrix.m());
2946 const number *matrix_ptr = &local_matrix(loc_row, 0);
2947
2948 for (size_type j = column_start; j < column_end; ++j)
2949 {
2950 const size_type loc_col = global_cols.local_row(j);
2951 AssertIndexRange(loc_col, local_matrix.n());
2952 const number col_val = matrix_ptr[loc_col];
2953 if (col_val != number())
2954 {
2955 *val_ptr++ = static_cast<number>(col_val);
2956 *col_ptr++ = global_cols.global_row(j);
2957 }
2958 }
2959 }
2960
2961 // more difficult part when there are indirect references and when we need
2962 // to do some more checks.
2963 else
2964 {
2965 for (size_type j = column_start; j < column_end; ++j)
2966 {
2967 number col_val = resolve_matrix_entry(
2968 global_rows, global_cols, i, j, loc_row, local_matrix);
2969
2970 // if we got some nontrivial value, append it to the array of
2971 // values.
2972 if (col_val != number())
2973 {
2974 *val_ptr++ = static_cast<number>(col_val);
2975 *col_ptr++ = global_cols.global_row(j);
2976 }
2977 }
2978 }
2979 }
2980
2981 // specialized function that can write into the row of a
2982 // SparseMatrix<number>.
2983 namespace dealiiSparseMatrix
2984 {
2985 template <typename SparseMatrixIterator, typename LocalType>
2986 static inline void
add_value(const LocalType value,const size_type row,const size_type column,SparseMatrixIterator & matrix_values)2987 add_value(const LocalType value,
2988 const size_type row,
2989 const size_type column,
2990 SparseMatrixIterator &matrix_values)
2991 {
2992 (void)row;
2993 if (value != LocalType())
2994 {
2995 while (matrix_values->column() < column)
2996 ++matrix_values;
2997 Assert(matrix_values->column() == column,
2998 typename SparseMatrix<
2999 typename SparseMatrixIterator::MatrixType::value_type>::
3000 ExcInvalidIndex(row, column));
3001 matrix_values->value() += value;
3002 }
3003 }
3004 } // namespace dealiiSparseMatrix
3005
3006 // similar as before, now with shortcut for deal.II sparse matrices. this
3007 // lets us avoid using extra arrays, and does all the operations just in
3008 // place, i.e., in the respective matrix row
3009 template <typename number>
3010 inline void
resolve_matrix_row(const GlobalRowsFromLocal<number> & global_rows,const size_type i,const size_type column_start,const size_type column_end,const FullMatrix<number> & local_matrix,SparseMatrix<number> * sparse_matrix)3011 resolve_matrix_row(const GlobalRowsFromLocal<number> &global_rows,
3012 const size_type i,
3013 const size_type column_start,
3014 const size_type column_end,
3015 const FullMatrix<number> & local_matrix,
3016 SparseMatrix<number> * sparse_matrix)
3017 {
3018 if (column_end == column_start)
3019 return;
3020
3021 AssertIndexRange(column_end - 1, global_rows.size());
3022 const SparsityPattern &sparsity = sparse_matrix->get_sparsity_pattern();
3023
3024 if (sparsity.n_nonzero_elements() == 0)
3025 return;
3026
3027 const size_type row = global_rows.global_row(i);
3028 const size_type loc_row = global_rows.local_row(i);
3029
3030 typename SparseMatrix<number>::iterator matrix_values =
3031 sparse_matrix->begin(row);
3032 const bool optimize_diagonal = sparsity.n_rows() == sparsity.n_cols();
3033
3034 // distinguish three cases about what can happen for checking whether the
3035 // diagonal is the first element of the row. this avoids if statements at
3036 // the innermost loop positions
3037
3038 if (!optimize_diagonal) // case 1: no diagonal optimization in matrix
3039 {
3040 if (global_rows.have_indirect_rows() == false)
3041 {
3042 AssertIndexRange(loc_row, local_matrix.m());
3043 const number *matrix_ptr = &local_matrix(loc_row, 0);
3044
3045 for (size_type j = column_start; j < column_end; ++j)
3046 {
3047 const size_type loc_col = global_rows.local_row(j);
3048 const number col_val = matrix_ptr[loc_col];
3049 dealiiSparseMatrix::add_value(col_val,
3050 row,
3051 global_rows.global_row(j),
3052 matrix_values);
3053 }
3054 }
3055 else
3056 {
3057 for (size_type j = column_start; j < column_end; ++j)
3058 {
3059 number col_val = resolve_matrix_entry(
3060 global_rows, global_rows, i, j, loc_row, local_matrix);
3061 dealiiSparseMatrix::add_value(col_val,
3062 row,
3063 global_rows.global_row(j),
3064 matrix_values);
3065 }
3066 }
3067 }
3068 else if (i >= column_start && i < column_end) // case 2: can split loop
3069 {
3070 ++matrix_values; // jump over diagonal element
3071 if (global_rows.have_indirect_rows() == false)
3072 {
3073 AssertIndexRange(loc_row, local_matrix.m());
3074 const number *matrix_ptr = &local_matrix(loc_row, 0);
3075
3076 sparse_matrix->begin(row)->value() += matrix_ptr[loc_row];
3077 for (size_type j = column_start; j < i; ++j)
3078 {
3079 const size_type loc_col = global_rows.local_row(j);
3080 const number col_val = matrix_ptr[loc_col];
3081 dealiiSparseMatrix::add_value(col_val,
3082 row,
3083 global_rows.global_row(j),
3084 matrix_values);
3085 }
3086 for (size_type j = i + 1; j < column_end; ++j)
3087 {
3088 const size_type loc_col = global_rows.local_row(j);
3089 const number col_val = matrix_ptr[loc_col];
3090 dealiiSparseMatrix::add_value(col_val,
3091 row,
3092 global_rows.global_row(j),
3093 matrix_values);
3094 }
3095 }
3096 else
3097 {
3098 sparse_matrix->begin(row)->value() += resolve_matrix_entry(
3099 global_rows, global_rows, i, i, loc_row, local_matrix);
3100 for (size_type j = column_start; j < i; ++j)
3101 {
3102 number col_val = resolve_matrix_entry(
3103 global_rows, global_rows, i, j, loc_row, local_matrix);
3104 dealiiSparseMatrix::add_value(col_val,
3105 row,
3106 global_rows.global_row(j),
3107 matrix_values);
3108 }
3109 for (size_type j = i + 1; j < column_end; ++j)
3110 {
3111 number col_val = resolve_matrix_entry(
3112 global_rows, global_rows, i, j, loc_row, local_matrix);
3113 dealiiSparseMatrix::add_value(col_val,
3114 row,
3115 global_rows.global_row(j),
3116 matrix_values);
3117 }
3118 }
3119 }
3120 // case 3: can't say - need to check inside the loop
3121 else if (global_rows.have_indirect_rows() == false)
3122 {
3123 ++matrix_values; // jump over diagonal element
3124 AssertIndexRange(loc_row, local_matrix.m());
3125 const number *matrix_ptr = &local_matrix(loc_row, 0);
3126
3127 for (size_type j = column_start; j < column_end; ++j)
3128 {
3129 const size_type loc_col = global_rows.local_row(j);
3130 const number col_val = matrix_ptr[loc_col];
3131 if (row == global_rows.global_row(j))
3132 sparse_matrix->begin(row)->value() += col_val;
3133 else
3134 dealiiSparseMatrix::add_value(col_val,
3135 row,
3136 global_rows.global_row(j),
3137 matrix_values);
3138 }
3139 }
3140 else
3141 {
3142 ++matrix_values; // jump over diagonal element
3143 for (size_type j = column_start; j < column_end; ++j)
3144 {
3145 number col_val = resolve_matrix_entry(
3146 global_rows, global_rows, i, j, loc_row, local_matrix);
3147 if (row == global_rows.global_row(j))
3148 sparse_matrix->begin(row)->value() += col_val;
3149 else
3150 dealiiSparseMatrix::add_value(col_val,
3151 row,
3152 global_rows.global_row(j),
3153 matrix_values);
3154 }
3155 }
3156 }
3157
3158 // Same function to resolve all entries that will be added to the given
3159 // global row global_rows[i] as before, now for sparsity pattern
3160 template <typename number>
3161 inline void
resolve_matrix_row(const GlobalRowsFromLocal<number> & global_rows,const size_type i,const size_type column_start,const size_type column_end,const Table<2,bool> & dof_mask,std::vector<size_type>::iterator & col_ptr)3162 resolve_matrix_row(const GlobalRowsFromLocal<number> &global_rows,
3163 const size_type i,
3164 const size_type column_start,
3165 const size_type column_end,
3166 const Table<2, bool> & dof_mask,
3167 std::vector<size_type>::iterator & col_ptr)
3168 {
3169 if (column_end == column_start)
3170 return;
3171
3172 const size_type loc_row = global_rows.local_row(i);
3173
3174 // fast function if there are no indirect references to any of the local
3175 // rows at all on this set of dofs
3176 if (global_rows.have_indirect_rows() == false)
3177 {
3178 Assert(loc_row < dof_mask.n_rows(), ExcInternalError());
3179
3180 for (size_type j = column_start; j < column_end; ++j)
3181 {
3182 const size_type loc_col = global_rows.local_row(j);
3183 Assert(loc_col < dof_mask.n_cols(), ExcInternalError());
3184
3185 if (dof_mask(loc_row, loc_col) == true)
3186 *col_ptr++ = global_rows.global_row(j);
3187 }
3188 }
3189
3190 // slower functions when there are indirect references and when we need to
3191 // do some more checks.
3192 else
3193 {
3194 for (size_type j = column_start; j < column_end; ++j)
3195 {
3196 const size_type loc_col = global_rows.local_row(j);
3197 if (loc_row != numbers::invalid_size_type)
3198 {
3199 Assert(loc_row < dof_mask.n_rows(), ExcInternalError());
3200 if (loc_col != numbers::invalid_size_type)
3201 {
3202 Assert(loc_col < dof_mask.n_cols(), ExcInternalError());
3203 if (dof_mask(loc_row, loc_col) == true)
3204 goto add_this_index;
3205 }
3206
3207 for (size_type p = 0; p < global_rows.size(j); ++p)
3208 if (dof_mask(loc_row, global_rows.local_row(j, p)) == true)
3209 goto add_this_index;
3210 }
3211
3212 for (size_type q = 0; q < global_rows.size(i); ++q)
3213 {
3214 if (loc_col != numbers::invalid_size_type)
3215 {
3216 Assert(loc_col < dof_mask.n_cols(), ExcInternalError());
3217 if (dof_mask(global_rows.local_row(i, q), loc_col) ==
3218 true)
3219 goto add_this_index;
3220 }
3221
3222 for (size_type p = 0; p < global_rows.size(j); ++p)
3223 if (dof_mask(global_rows.local_row(i, q),
3224 global_rows.local_row(j, p)) == true)
3225 goto add_this_index;
3226 }
3227
3228 continue;
3229 // if we got some nontrivial value, append it to the array of
3230 // values.
3231 add_this_index:
3232 *col_ptr++ = global_rows.global_row(j);
3233 }
3234 }
3235 }
3236
3237 // to make sure that the global matrix remains invertible, we need to do
3238 // something with the diagonal elements. Add the average of the
3239 // absolute values of the local matrix diagonals, so the resulting entry
3240 // will always be positive and furthermore be in the same order of magnitude
3241 // as the other elements of the matrix. If all local matrix diagonals are
3242 // zero, add the l1 norm of the local matrix divided by the matrix size
3243 // to the diagonal of the global matrix. If the entire local matrix is zero,
3244 // add 1 to the diagonal of the global matrix.
3245 //
3246 // note that this also captures the special case that a dof is both
3247 // constrained and fixed (this can happen for hanging nodes in 3d that also
3248 // happen to be on the boundary). in that case, following the program flow
3249 // in distribute_local_to_global, it is realized that when distributing the
3250 // row and column no elements of the matrix are actually touched if all the
3251 // degrees of freedom to which this dof is constrained are also constrained
3252 // (the usual case with hanging nodes in 3d). however, in the line below, we
3253 // do actually do something with this dof
3254 template <typename number, typename MatrixType, typename VectorType>
3255 inline void
set_matrix_diagonals(const internal::AffineConstraints::GlobalRowsFromLocal<number> & global_rows,const std::vector<size_type> & local_dof_indices,const FullMatrix<number> & local_matrix,const dealii::AffineConstraints<number> & constraints,MatrixType & global_matrix,VectorType & global_vector,bool use_inhomogeneities_for_rhs)3256 set_matrix_diagonals(
3257 const internal::AffineConstraints::GlobalRowsFromLocal<number>
3258 & global_rows,
3259 const std::vector<size_type> & local_dof_indices,
3260 const FullMatrix<number> & local_matrix,
3261 const dealii::AffineConstraints<number> &constraints,
3262 MatrixType & global_matrix,
3263 VectorType & global_vector,
3264 bool use_inhomogeneities_for_rhs)
3265 {
3266 if (global_rows.n_constraints() > 0)
3267 {
3268 number average_diagonal = number();
3269 for (size_type i = 0; i < local_matrix.m(); ++i)
3270 average_diagonal += std::abs(local_matrix(i, i));
3271 average_diagonal /= static_cast<number>(local_matrix.m());
3272
3273 // handle the case that all diagonal elements are zero
3274 if (average_diagonal == static_cast<number>(0.))
3275 {
3276 average_diagonal = static_cast<number>(local_matrix.l1_norm()) /
3277 static_cast<number>(local_matrix.m());
3278 // if the entire matrix is zero, use 1. for the diagonal
3279 if (average_diagonal == static_cast<number>(0.))
3280 average_diagonal = static_cast<number>(1.);
3281 }
3282
3283 for (size_type i = 0; i < global_rows.n_constraints(); i++)
3284 {
3285 const size_type local_row = global_rows.constraint_origin(i);
3286 const size_type global_row = local_dof_indices[local_row];
3287
3288 const number current_diagonal =
3289 local_matrix(local_row, local_row);
3290 if (std::abs(current_diagonal) != 0.)
3291 {
3292 global_matrix.add(global_row,
3293 global_row,
3294 std::abs(current_diagonal));
3295 // if the use_inhomogeneities_for_rhs flag is set to true, the
3296 // inhomogeneities are used to create the global vector.
3297 // instead of fill in a zero in the ith components with an
3298 // inhomogeneity, we set those to:
3299 // inhomogeneity(i)*global_matrix (i,i).
3300 if (use_inhomogeneities_for_rhs == true)
3301 global_vector(global_row) +=
3302 current_diagonal *
3303 constraints.get_inhomogeneity(global_row);
3304 }
3305 else
3306 {
3307 global_matrix.add(global_row, global_row, average_diagonal);
3308
3309 if (use_inhomogeneities_for_rhs == true)
3310 global_vector(global_row) +=
3311 average_diagonal *
3312 constraints.get_inhomogeneity(global_row);
3313 }
3314 }
3315 }
3316 }
3317
3318 // similar function as the one above for setting matrix diagonals, but now
3319 // doing that for sparsity patterns when setting them up using
3320 // add_entries_local_to_global. In case we keep constrained entries, add all
3321 // the rows and columns related to the constrained dof, otherwise just add
3322 // the diagonal
3323 template <typename number, typename SparsityPatternType>
3324 inline void
set_sparsity_diagonals(const internal::AffineConstraints::GlobalRowsFromLocal<number> & global_rows,const std::vector<size_type> & local_dof_indices,const Table<2,bool> & dof_mask,const bool keep_constrained_entries,SparsityPatternType & sparsity_pattern)3325 set_sparsity_diagonals(
3326 const internal::AffineConstraints::GlobalRowsFromLocal<number>
3327 & global_rows,
3328 const std::vector<size_type> &local_dof_indices,
3329 const Table<2, bool> & dof_mask,
3330 const bool keep_constrained_entries,
3331 SparsityPatternType & sparsity_pattern)
3332 {
3333 // if we got constraints, need to add the diagonal element and, if the
3334 // user requested so, also the rest of the entries in rows and columns
3335 // that have been left out above
3336 if (global_rows.n_constraints() > 0)
3337 {
3338 for (size_type i = 0; i < global_rows.n_constraints(); i++)
3339 {
3340 const size_type local_row = global_rows.constraint_origin(i);
3341 const size_type global_row = local_dof_indices[local_row];
3342 if (keep_constrained_entries == true)
3343 {
3344 for (size_type j = 0; j < local_dof_indices.size(); ++j)
3345 {
3346 if (dof_mask(local_row, j) == true)
3347 sparsity_pattern.add(global_row, local_dof_indices[j]);
3348 if (dof_mask(j, local_row) == true)
3349 sparsity_pattern.add(local_dof_indices[j], global_row);
3350 }
3351 }
3352 else
3353 // don't keep constrained entries - just add the diagonal.
3354 sparsity_pattern.add(global_row, global_row);
3355 }
3356 }
3357 }
3358
3359 } // end of namespace AffineConstraints
3360 } // end of namespace internal
3361
3362 // Basic idea of setting up a list of
3363 // all global dofs: first find all rows and columns
3364 // that we are going to write touch,
3365 // and then go through the
3366 // lines and collect all the local rows that
3367 // are related to it.
3368 template <typename number>
3369 void
make_sorted_row_list(const std::vector<size_type> & local_dof_indices,internal::AffineConstraints::GlobalRowsFromLocal<number> & global_rows)3370 AffineConstraints<number>::make_sorted_row_list(
3371 const std::vector<size_type> & local_dof_indices,
3372 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows) const
3373 {
3374 const size_type n_local_dofs = local_dof_indices.size();
3375 AssertDimension(n_local_dofs, global_rows.size());
3376
3377 // when distributing the local data to the global matrix, we can quite
3378 // cheaply sort the indices (obviously, this introduces the need for
3379 // allocating some memory on the way, but we need to do this only for rows,
3380 // whereas the distribution process itself goes over rows and columns). This
3381 // has the advantage that when writing into the global matrix, we can make
3382 // use of the sortedness.
3383
3384 // so the first step is to create a sorted list of all row values that are
3385 // possible. these values are either the rows from unconstrained dofs, or
3386 // some indices introduced by dofs constrained to a combination of some
3387 // other dofs. regarding the data type, choose a <tt>std::vector</tt> of a
3388 // pair of unsigned ints (for global columns) and internal data (containing
3389 // local columns + possible jumps from constraints). Choosing
3390 // <tt>std::map</tt> or anything else M.K. knows of would be much more
3391 // expensive here!
3392
3393 // cache whether we have to resolve any indirect rows generated from
3394 // resolving constrained dofs.
3395 size_type added_rows = 0;
3396
3397 // first add the indices in an unsorted way and only keep track of the
3398 // constraints that appear. They are resolved in a second step.
3399 for (size_type i = 0; i < n_local_dofs; ++i)
3400 {
3401 if (is_constrained(local_dof_indices[i]) == false)
3402 {
3403 global_rows.global_row(added_rows) = local_dof_indices[i];
3404 global_rows.local_row(added_rows++) = i;
3405 }
3406 else
3407 global_rows.insert_constraint(i);
3408 }
3409 global_rows.sort();
3410
3411 const size_type n_constrained_rows = n_local_dofs - added_rows;
3412 for (size_type i = 0; i < n_constrained_rows; ++i)
3413 {
3414 const size_type local_row = global_rows.constraint_origin(i);
3415 AssertIndexRange(local_row, n_local_dofs);
3416 const size_type global_row = local_dof_indices[local_row];
3417 Assert(is_constrained(global_row), ExcInternalError());
3418 const ConstraintLine &position =
3419 lines[lines_cache[calculate_line_index(global_row)]];
3420 if (position.inhomogeneity != number(0.))
3421 global_rows.set_ith_constraint_inhomogeneous(i);
3422 for (size_type q = 0; q < position.entries.size(); ++q)
3423 global_rows.insert_index(position.entries[q].first,
3424 local_row,
3425 position.entries[q].second);
3426 }
3427 }
3428
3429 // Same function as before, but now do only extract the global indices that
3430 // come from the local ones without storing their origin. Used for sparsity
3431 // pattern generation.
3432 template <typename number>
3433 inline void
make_sorted_row_list(const std::vector<size_type> & local_dof_indices,std::vector<size_type> & active_dofs)3434 AffineConstraints<number>::make_sorted_row_list(
3435 const std::vector<size_type> &local_dof_indices,
3436 std::vector<size_type> & active_dofs) const
3437 {
3438 const size_type n_local_dofs = local_dof_indices.size();
3439 size_type added_rows = 0;
3440 for (size_type i = 0; i < n_local_dofs; ++i)
3441 {
3442 if (is_constrained(local_dof_indices[i]) == false)
3443 {
3444 active_dofs[added_rows++] = local_dof_indices[i];
3445 continue;
3446 }
3447
3448 active_dofs[n_local_dofs - i + added_rows - 1] = i;
3449 }
3450 std::sort(active_dofs.begin(), active_dofs.begin() + added_rows);
3451
3452 const size_type n_constrained_dofs = n_local_dofs - added_rows;
3453 for (size_type i = n_constrained_dofs; i > 0; --i)
3454 {
3455 const size_type local_row = active_dofs.back();
3456
3457 // remove constrained entry since we are going to resolve it in place
3458 active_dofs.pop_back();
3459 const size_type global_row = local_dof_indices[local_row];
3460 const ConstraintLine &position =
3461 lines[lines_cache[calculate_line_index(global_row)]];
3462 for (size_type q = 0; q < position.entries.size(); ++q)
3463 {
3464 const size_type new_index = position.entries[q].first;
3465 if (active_dofs[active_dofs.size() - i] < new_index)
3466 active_dofs.insert(active_dofs.end() - i + 1, new_index);
3467
3468 // make binary search to find where to put the new index in order to
3469 // keep the list sorted
3470 else
3471 {
3472 std::vector<size_type>::iterator it =
3473 Utilities::lower_bound(active_dofs.begin(),
3474 active_dofs.end() - i + 1,
3475 new_index);
3476 if (*it != new_index)
3477 active_dofs.insert(it, new_index);
3478 }
3479 }
3480 }
3481 }
3482
3483 // Resolve the constraints from the vector and apply inhomogeneities.
3484 template <typename number>
3485 template <typename MatrixScalar, typename VectorScalar>
3486 inline typename ProductType<VectorScalar, MatrixScalar>::type
resolve_vector_entry(const size_type i,const internal::AffineConstraints::GlobalRowsFromLocal<number> & global_rows,const Vector<VectorScalar> & local_vector,const std::vector<size_type> & local_dof_indices,const FullMatrix<MatrixScalar> & local_matrix)3487 AffineConstraints<number>::resolve_vector_entry(
3488 const size_type i,
3489 const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
3490 const Vector<VectorScalar> & local_vector,
3491 const std::vector<size_type> & local_dof_indices,
3492 const FullMatrix<MatrixScalar> &local_matrix) const
3493 {
3494 const size_type loc_row = global_rows.local_row(i);
3495 const size_type n_inhomogeneous_rows = global_rows.n_inhomogeneities();
3496 typename ProductType<VectorScalar, MatrixScalar>::type val = 0;
3497 // has a direct contribution from some local entry. If we have inhomogeneous
3498 // constraints, compute the contribution of the inhomogeneity in the current
3499 // row.
3500 if (loc_row != numbers::invalid_size_type)
3501 {
3502 val = local_vector(loc_row);
3503 for (size_type i = 0; i < n_inhomogeneous_rows; ++i)
3504 val -= (local_matrix(loc_row, global_rows.constraint_origin(i)) *
3505 lines[lines_cache[calculate_line_index(
3506 local_dof_indices[global_rows.constraint_origin(i)])]]
3507 .inhomogeneity);
3508 }
3509
3510 // go through the indirect contributions
3511 for (size_type q = 0; q < global_rows.size(i); ++q)
3512 {
3513 const size_type loc_row_q = global_rows.local_row(i, q);
3514 typename ProductType<VectorScalar, MatrixScalar>::type add_this =
3515 local_vector(loc_row_q);
3516 for (size_type k = 0; k < n_inhomogeneous_rows; ++k)
3517 add_this -=
3518 (local_matrix(loc_row_q, global_rows.constraint_origin(k)) *
3519 lines[lines_cache[calculate_line_index(
3520 local_dof_indices[global_rows.constraint_origin(k)])]]
3521 .inhomogeneity);
3522 val += add_this * global_rows.constraint_value(i, q);
3523 }
3524 return val;
3525 }
3526
3527 // internal implementation for distribute_local_to_global for standard
3528 // (non-block) matrices
3529 template <typename number>
3530 template <typename MatrixType, typename VectorType>
3531 void
distribute_local_to_global(const FullMatrix<number> & local_matrix,const Vector<number> & local_vector,const std::vector<size_type> & local_dof_indices,MatrixType & global_matrix,VectorType & global_vector,bool use_inhomogeneities_for_rhs,std::integral_constant<bool,false>)3532 AffineConstraints<number>::distribute_local_to_global(
3533 const FullMatrix<number> & local_matrix,
3534 const Vector<number> & local_vector,
3535 const std::vector<size_type> &local_dof_indices,
3536 MatrixType & global_matrix,
3537 VectorType & global_vector,
3538 bool use_inhomogeneities_for_rhs,
3539 std::integral_constant<bool, false>) const
3540 {
3541 // FIXME: static_assert MatrixType::value_type == number
3542
3543 // check whether we work on real vectors or we just used a dummy when
3544 // calling the other function above.
3545 const bool use_vectors =
3546 (local_vector.size() == 0 && global_vector.size() == 0) ? false : true;
3547 const bool use_dealii_matrix =
3548 std::is_same<MatrixType, SparseMatrix<number>>::value;
3549
3550 AssertDimension(local_matrix.n(), local_dof_indices.size());
3551 AssertDimension(local_matrix.m(), local_dof_indices.size());
3552 Assert(global_matrix.m() == global_matrix.n(), ExcNotQuadratic());
3553 if (use_vectors == true)
3554 {
3555 AssertDimension(local_matrix.m(), local_vector.size());
3556 AssertDimension(global_matrix.m(), global_vector.size());
3557 }
3558 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
3559
3560 const size_type n_local_dofs = local_dof_indices.size();
3561
3562 typename internal::AffineConstraints::ScratchDataAccessor<number>
3563 scratch_data(this->scratch_data);
3564
3565 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows =
3566 scratch_data->global_rows;
3567 global_rows.reinit(n_local_dofs);
3568 make_sorted_row_list(local_dof_indices, global_rows);
3569
3570 const size_type n_actual_dofs = global_rows.size();
3571
3572 // create arrays for the column data (indices and values) that will then be
3573 // written into the matrix. Shortcut for deal.II sparse matrix. We can use
3574 // the scratch data if we have a double matrix. Otherwise, we need to create
3575 // an array in any case since we cannot know about the actual data type in
3576 // the AffineConstraints class (unless we do cast). This involves a little
3577 // bit of logic to determine the type of the matrix value.
3578 std::vector<size_type> &cols = scratch_data->columns;
3579 std::vector<number> & vals = scratch_data->values;
3580 // create arrays for writing into the vector as well
3581 std::vector<size_type> &vector_indices = scratch_data->vector_indices;
3582 std::vector<typename VectorType::value_type> &vector_values =
3583 scratch_data->vector_values;
3584 vector_indices.resize(n_actual_dofs);
3585 vector_values.resize(n_actual_dofs);
3586 SparseMatrix<number> *sparse_matrix =
3587 dynamic_cast<SparseMatrix<number> *>(&global_matrix);
3588 if (use_dealii_matrix == false)
3589 {
3590 cols.resize(n_actual_dofs);
3591 vals.resize(n_actual_dofs);
3592 }
3593 else
3594 Assert(sparse_matrix != nullptr, ExcInternalError());
3595
3596 // now do the actual job. go through all the global rows that we will touch
3597 // and call resolve_matrix_row for each of those.
3598 size_type local_row_n = 0;
3599 for (size_type i = 0; i < n_actual_dofs; ++i)
3600 {
3601 const size_type row = global_rows.global_row(i);
3602
3603 // calculate all the data that will be written into the matrix row.
3604 if (use_dealii_matrix == false)
3605 {
3606 size_type *col_ptr = cols.data();
3607 // cast is uncritical here and only used to avoid compiler
3608 // warnings. We never access a non-double array
3609 number *val_ptr = vals.data();
3610 internal::AffineConstraints::resolve_matrix_row(global_rows,
3611 global_rows,
3612 i,
3613 0,
3614 n_actual_dofs,
3615 local_matrix,
3616 col_ptr,
3617 val_ptr);
3618 const size_type n_values = col_ptr - cols.data();
3619 if (n_values > 0)
3620 global_matrix.add(
3621 row, n_values, cols.data(), vals.data(), false, true);
3622 }
3623 else
3624 internal::AffineConstraints::resolve_matrix_row(
3625 global_rows, i, 0, n_actual_dofs, local_matrix, sparse_matrix);
3626
3627 // now to the vectors. besides doing the same job as we did above (i.e.,
3628 // distribute the content of the local vector into the global one), need
3629 // to account for inhomogeneities here: this corresponds to eliminating
3630 // the respective column in the local matrix with value on the right
3631 // hand side.
3632 if (use_vectors == true)
3633 {
3634 const typename VectorType::value_type val = resolve_vector_entry(
3635 i, global_rows, local_vector, local_dof_indices, local_matrix);
3636 AssertIsFinite(val);
3637
3638 if (val != typename VectorType::value_type())
3639 {
3640 vector_indices[local_row_n] = row;
3641 vector_values[local_row_n] = val;
3642 ++local_row_n;
3643 }
3644 }
3645 }
3646 // Drop the elements of vector_indices and vector_values that we do not use
3647 // (we may always elide writing zero values to vectors)
3648 const size_type n_local_rows = local_row_n;
3649 vector_indices.resize(n_local_rows);
3650 vector_values.resize(n_local_rows);
3651
3652 // While the standard case is that these types are equal, they need not be, so
3653 // only do a bulk update if they are. Note that the types in the arguments to
3654 // add must be equal if we have a Trilinos or PETSc vector but do not have to
3655 // be if we have a deal.II native vector: one could further optimize this for
3656 // Vector, LinearAlgebra::distributed::vector, etc.
3657 if (std::is_same<typename VectorType::value_type, number>::value)
3658 {
3659 global_vector.add(vector_indices,
3660 *reinterpret_cast<std::vector<number> *>(
3661 &vector_values));
3662 }
3663 else
3664 {
3665 for (size_type row_n = 0; row_n < n_local_rows; ++row_n)
3666 {
3667 global_vector(vector_indices[row_n]) +=
3668 static_cast<typename VectorType::value_type>(vector_values[row_n]);
3669 }
3670 }
3671
3672 internal::AffineConstraints::set_matrix_diagonals(
3673 global_rows,
3674 local_dof_indices,
3675 local_matrix,
3676 *this,
3677 global_matrix,
3678 global_vector,
3679 use_inhomogeneities_for_rhs);
3680 }
3681
3682 // similar function as above, but now specialized for block matrices. See the
3683 // other function for additional comments.
3684 template <typename number>
3685 template <typename MatrixType, typename VectorType>
3686 void
distribute_local_to_global(const FullMatrix<number> & local_matrix,const Vector<number> & local_vector,const std::vector<size_type> & local_dof_indices,MatrixType & global_matrix,VectorType & global_vector,bool use_inhomogeneities_for_rhs,std::integral_constant<bool,true>)3687 AffineConstraints<number>::distribute_local_to_global(
3688 const FullMatrix<number> & local_matrix,
3689 const Vector<number> & local_vector,
3690 const std::vector<size_type> &local_dof_indices,
3691 MatrixType & global_matrix,
3692 VectorType & global_vector,
3693 bool use_inhomogeneities_for_rhs,
3694 std::integral_constant<bool, true>) const
3695 {
3696 const bool use_vectors =
3697 (local_vector.size() == 0 && global_vector.size() == 0) ? false : true;
3698 const bool use_dealii_matrix =
3699 std::is_same<MatrixType, BlockSparseMatrix<number>>::value;
3700
3701 AssertDimension(local_matrix.n(), local_dof_indices.size());
3702 AssertDimension(local_matrix.m(), local_dof_indices.size());
3703 Assert(global_matrix.m() == global_matrix.n(), ExcNotQuadratic());
3704 Assert(global_matrix.n_block_rows() == global_matrix.n_block_cols(),
3705 ExcNotQuadratic());
3706 if (use_vectors == true)
3707 {
3708 AssertDimension(local_matrix.m(), local_vector.size());
3709 AssertDimension(global_matrix.m(), global_vector.size());
3710 }
3711 Assert(sorted == true, ExcMatrixNotClosed());
3712
3713 typename internal::AffineConstraints::ScratchDataAccessor<number>
3714 scratch_data(this->scratch_data);
3715
3716 const size_type n_local_dofs = local_dof_indices.size();
3717 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows =
3718 scratch_data->global_rows;
3719 global_rows.reinit(n_local_dofs);
3720
3721 make_sorted_row_list(local_dof_indices, global_rows);
3722 const size_type n_actual_dofs = global_rows.size();
3723
3724 std::vector<size_type> &global_indices = scratch_data->vector_indices;
3725 if (use_vectors == true)
3726 {
3727 global_indices.resize(n_actual_dofs);
3728 for (size_type i = 0; i < n_actual_dofs; ++i)
3729 global_indices[i] = global_rows.global_row(i);
3730 }
3731
3732 // additional construct that also takes care of block indices.
3733 const size_type num_blocks = global_matrix.n_block_rows();
3734 std::vector<size_type> &block_starts = scratch_data->block_starts;
3735 block_starts.resize(num_blocks + 1);
3736 internal::AffineConstraints::make_block_starts(global_matrix,
3737 global_rows,
3738 block_starts);
3739
3740 std::vector<size_type> &cols = scratch_data->columns;
3741 std::vector<number> & vals = scratch_data->values;
3742 if (use_dealii_matrix == false)
3743 {
3744 cols.resize(n_actual_dofs);
3745 vals.resize(n_actual_dofs);
3746 }
3747
3748 // the basic difference to the non-block variant from now onwards is that we
3749 // go through the blocks of the matrix separately, which allows us to set
3750 // the block entries individually
3751 for (size_type block = 0; block < num_blocks; ++block)
3752 {
3753 const size_type next_block = block_starts[block + 1];
3754 for (size_type i = block_starts[block]; i < next_block; ++i)
3755 {
3756 const size_type row = global_rows.global_row(i);
3757
3758 for (size_type block_col = 0; block_col < num_blocks; ++block_col)
3759 {
3760 const size_type start_block = block_starts[block_col],
3761 end_block = block_starts[block_col + 1];
3762 if (use_dealii_matrix == false)
3763 {
3764 size_type *col_ptr = cols.data();
3765 number * val_ptr = vals.data();
3766 internal::AffineConstraints::resolve_matrix_row(global_rows,
3767 global_rows,
3768 i,
3769 start_block,
3770 end_block,
3771 local_matrix,
3772 col_ptr,
3773 val_ptr);
3774 const size_type n_values = col_ptr - cols.data();
3775 if (n_values > 0)
3776 global_matrix.block(block, block_col)
3777 .add(
3778 row, n_values, cols.data(), vals.data(), false, true);
3779 }
3780 else
3781 {
3782 SparseMatrix<number> *sparse_matrix =
3783 dynamic_cast<SparseMatrix<number> *>(
3784 &global_matrix.block(block, block_col));
3785 Assert(sparse_matrix != nullptr, ExcInternalError());
3786 internal::AffineConstraints::resolve_matrix_row(
3787 global_rows,
3788 i,
3789 start_block,
3790 end_block,
3791 local_matrix,
3792 sparse_matrix);
3793 }
3794 }
3795
3796 if (use_vectors == true)
3797 {
3798 const number val = resolve_vector_entry(
3799 i, global_rows, local_vector, local_dof_indices, local_matrix);
3800
3801 if (val != number())
3802 global_vector(global_indices[i]) +=
3803 static_cast<typename VectorType::value_type>(val);
3804 }
3805 }
3806 }
3807
3808 internal::AffineConstraints::set_matrix_diagonals(
3809 global_rows,
3810 local_dof_indices,
3811 local_matrix,
3812 *this,
3813 global_matrix,
3814 global_vector,
3815 use_inhomogeneities_for_rhs);
3816 }
3817
3818
3819 template <typename number>
3820 template <typename MatrixType>
3821 void
distribute_local_to_global(const FullMatrix<number> & local_matrix,const std::vector<size_type> & row_indices,const std::vector<size_type> & col_indices,MatrixType & global_matrix)3822 AffineConstraints<number>::distribute_local_to_global(
3823 const FullMatrix<number> & local_matrix,
3824 const std::vector<size_type> &row_indices,
3825 const std::vector<size_type> &col_indices,
3826 MatrixType & global_matrix) const
3827 {
3828 distribute_local_to_global(
3829 local_matrix, row_indices, *this, col_indices, global_matrix);
3830 }
3831
3832
3833 template <typename number>
3834 template <typename MatrixType>
3835 void
distribute_local_to_global(const FullMatrix<number> & local_matrix,const std::vector<size_type> & row_indices,const AffineConstraints<number> & col_constraint_matrix,const std::vector<size_type> & col_indices,MatrixType & global_matrix)3836 AffineConstraints<number>::distribute_local_to_global(
3837 const FullMatrix<number> & local_matrix,
3838 const std::vector<size_type> & row_indices,
3839 const AffineConstraints<number> &col_constraint_matrix,
3840 const std::vector<size_type> & col_indices,
3841 MatrixType & global_matrix) const
3842 {
3843 AssertDimension(local_matrix.m(), row_indices.size());
3844 AssertDimension(local_matrix.n(), col_indices.size());
3845
3846 const size_type n_local_row_dofs = row_indices.size();
3847 const size_type n_local_col_dofs = col_indices.size();
3848
3849 typename internal::AffineConstraints::ScratchDataAccessor<
3850 typename MatrixType::value_type>
3851 scratch_data(this->scratch_data);
3852
3853 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows =
3854 scratch_data->global_rows;
3855 global_rows.reinit(n_local_row_dofs);
3856
3857 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_cols =
3858 scratch_data->global_columns;
3859 global_cols.reinit(n_local_col_dofs);
3860
3861 make_sorted_row_list(row_indices, global_rows);
3862 col_constraint_matrix.make_sorted_row_list(col_indices, global_cols);
3863
3864 const size_type n_actual_row_dofs = global_rows.size();
3865 const size_type n_actual_col_dofs = global_cols.size();
3866
3867 // create arrays for the column data (indices and values) that will then be
3868 // written into the matrix. Shortcut for deal.II sparse matrix
3869 std::vector<size_type> &cols = scratch_data->columns;
3870 std::vector<number> & vals = scratch_data->values;
3871 cols.resize(n_actual_col_dofs);
3872 vals.resize(n_actual_col_dofs);
3873
3874 // now do the actual job.
3875 for (size_type i = 0; i < n_actual_row_dofs; ++i)
3876 {
3877 const size_type row = global_rows.global_row(i);
3878
3879 // calculate all the data that will be written into the matrix row.
3880 size_type *col_ptr = cols.data();
3881 number * val_ptr = vals.data();
3882 internal::AffineConstraints::resolve_matrix_row(global_rows,
3883 global_cols,
3884 i,
3885 0,
3886 n_actual_col_dofs,
3887 local_matrix,
3888 col_ptr,
3889 val_ptr);
3890 const size_type n_values = col_ptr - cols.data();
3891 if (n_values > 0)
3892 global_matrix.add(row, n_values, cols.data(), vals.data(), false, true);
3893 }
3894 }
3895
3896 template <typename number>
3897 template <typename SparsityPatternType>
3898 void
add_entries_local_to_global(const std::vector<size_type> & local_dof_indices,SparsityPatternType & sparsity_pattern,const bool keep_constrained_entries,const Table<2,bool> & dof_mask,std::integral_constant<bool,false>)3899 AffineConstraints<number>::add_entries_local_to_global(
3900 const std::vector<size_type> &local_dof_indices,
3901 SparsityPatternType & sparsity_pattern,
3902 const bool keep_constrained_entries,
3903 const Table<2, bool> & dof_mask,
3904 std::integral_constant<bool, false>) const
3905 {
3906 Assert(sparsity_pattern.n_rows() == sparsity_pattern.n_cols(),
3907 ExcNotQuadratic());
3908
3909 const size_type n_local_dofs = local_dof_indices.size();
3910 typename internal::AffineConstraints::ScratchDataAccessor<number>
3911 scratch_data(this->scratch_data);
3912
3913 const bool dof_mask_is_active = (dof_mask.n_rows() == n_local_dofs);
3914 if (dof_mask_is_active == true)
3915 {
3916 AssertDimension(dof_mask.n_cols(), n_local_dofs);
3917 }
3918 else
3919 {
3920 // if the dof mask is not active, all we have to do is to add some indices
3921 // in a matrix format. To do this, we first create an array of all the
3922 // indices that are to be added. these indices are the local dof indices
3923 // plus some indices that come from constraints.
3924 std::vector<size_type> &actual_dof_indices = scratch_data->columns;
3925 actual_dof_indices.resize(n_local_dofs);
3926 make_sorted_row_list(local_dof_indices, actual_dof_indices);
3927 const size_type n_actual_dofs = actual_dof_indices.size();
3928
3929 // now add the indices we collected above to the sparsity pattern. Very
3930 // easy here - just add the same array to all the rows...
3931 for (size_type i = 0; i < n_actual_dofs; ++i)
3932 sparsity_pattern.add_entries(actual_dof_indices[i],
3933 actual_dof_indices.begin(),
3934 actual_dof_indices.end(),
3935 true);
3936
3937 // need to add the whole row and column structure in case we keep
3938 // constrained entries. Unfortunately, we can't use the nice matrix
3939 // structure we use elsewhere, so manually add those indices one by one.
3940 for (size_type i = 0; i < n_local_dofs; i++)
3941 if (is_constrained(local_dof_indices[i]))
3942 {
3943 if (keep_constrained_entries == true)
3944 for (size_type j = 0; j < n_local_dofs; j++)
3945 {
3946 sparsity_pattern.add(local_dof_indices[i],
3947 local_dof_indices[j]);
3948 sparsity_pattern.add(local_dof_indices[j],
3949 local_dof_indices[i]);
3950 }
3951 else
3952 sparsity_pattern.add(local_dof_indices[i], local_dof_indices[i]);
3953 }
3954
3955 return;
3956 }
3957
3958 // complicated case: we need to filter out some indices. then the function
3959 // gets similar to the function for distributing matrix entries, see there
3960 // for additional comments.
3961 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows =
3962 scratch_data->global_rows;
3963 global_rows.reinit(n_local_dofs);
3964 make_sorted_row_list(local_dof_indices, global_rows);
3965 const size_type n_actual_dofs = global_rows.size();
3966
3967 // create arrays for the column indices that will then be written into the
3968 // sparsity pattern.
3969 std::vector<size_type> &cols = scratch_data->columns;
3970 cols.resize(n_actual_dofs);
3971
3972 for (size_type i = 0; i < n_actual_dofs; ++i)
3973 {
3974 std::vector<size_type>::iterator col_ptr = cols.begin();
3975 const size_type row = global_rows.global_row(i);
3976 internal::AffineConstraints::resolve_matrix_row(
3977 global_rows, i, 0, n_actual_dofs, dof_mask, col_ptr);
3978
3979 // finally, write all the information that accumulated under the given
3980 // process into the global matrix row and into the vector
3981 if (col_ptr != cols.begin())
3982 sparsity_pattern.add_entries(row, cols.begin(), col_ptr, true);
3983 }
3984 internal::AffineConstraints::set_sparsity_diagonals(global_rows,
3985 local_dof_indices,
3986 dof_mask,
3987 keep_constrained_entries,
3988 sparsity_pattern);
3989 }
3990
3991 template <typename number>
3992 template <typename SparsityPatternType>
3993 void
add_entries_local_to_global(const std::vector<size_type> & row_indices,const std::vector<size_type> & col_indices,SparsityPatternType & sparsity_pattern,const bool keep_constrained_entries,const Table<2,bool> & dof_mask)3994 AffineConstraints<number>::add_entries_local_to_global(
3995 const std::vector<size_type> &row_indices,
3996 const std::vector<size_type> &col_indices,
3997 SparsityPatternType & sparsity_pattern,
3998 const bool keep_constrained_entries,
3999 const Table<2, bool> & dof_mask) const
4000 {
4001 const size_type n_local_rows = row_indices.size();
4002 const size_type n_local_cols = col_indices.size();
4003
4004 // if constrained entries should be kept, need to add rows and columns of
4005 // those to the sparsity pattern
4006 if (keep_constrained_entries == true)
4007 {
4008 for (const size_type row_index : row_indices)
4009 if (is_constrained(row_index))
4010 for (const size_type col_index : col_indices)
4011 sparsity_pattern.add(row_index, col_index);
4012 for (const size_type col_index : col_indices)
4013 if (is_constrained(col_index))
4014 for (const size_type row_index : row_indices)
4015 sparsity_pattern.add(row_index, col_index);
4016 }
4017
4018 // if the dof mask is not active, all we have to do is to add some indices
4019 // in a matrix format. To do this, we first create an array of all the
4020 // indices that are to be added. these indices are the local dof indices
4021 // plus some indices that come from constraints.
4022 const bool dof_mask_is_active =
4023 dof_mask.n_rows() == n_local_rows && dof_mask.n_cols() == n_local_cols;
4024 if (dof_mask_is_active == false)
4025 {
4026 std::vector<size_type> actual_row_indices(n_local_rows);
4027 std::vector<size_type> actual_col_indices(n_local_cols);
4028 make_sorted_row_list(row_indices, actual_row_indices);
4029 make_sorted_row_list(col_indices, actual_col_indices);
4030 const size_type n_actual_rows = actual_row_indices.size();
4031
4032 // now add the indices we collected above to the sparsity pattern. Very
4033 // easy here - just add the same array to all the rows...
4034 for (size_type i = 0; i < n_actual_rows; ++i)
4035 sparsity_pattern.add_entries(actual_row_indices[i],
4036 actual_col_indices.begin(),
4037 actual_col_indices.end(),
4038 true);
4039 return;
4040 }
4041
4042 // TODO: implement this
4043 Assert(false, ExcNotImplemented());
4044 }
4045
4046 template <typename number>
4047 template <typename SparsityPatternType>
4048 void
add_entries_local_to_global(const std::vector<size_type> & local_dof_indices,SparsityPatternType & sparsity_pattern,const bool keep_constrained_entries,const Table<2,bool> & dof_mask,std::integral_constant<bool,true>)4049 AffineConstraints<number>::add_entries_local_to_global(
4050 const std::vector<size_type> &local_dof_indices,
4051 SparsityPatternType & sparsity_pattern,
4052 const bool keep_constrained_entries,
4053 const Table<2, bool> & dof_mask,
4054 std::integral_constant<bool, true>) const
4055 {
4056 // just as the other add_entries_local_to_global function, but now
4057 // specialized for block matrices.
4058 Assert(sparsity_pattern.n_rows() == sparsity_pattern.n_cols(),
4059 ExcNotQuadratic());
4060 Assert(sparsity_pattern.n_block_rows() == sparsity_pattern.n_block_cols(),
4061 ExcNotQuadratic());
4062
4063 const size_type n_local_dofs = local_dof_indices.size();
4064 const size_type num_blocks = sparsity_pattern.n_block_rows();
4065
4066 typename internal::AffineConstraints::ScratchDataAccessor<number>
4067 scratch_data(this->scratch_data);
4068
4069 const bool dof_mask_is_active = (dof_mask.n_rows() == n_local_dofs);
4070 if (dof_mask_is_active == true)
4071 {
4072 AssertDimension(dof_mask.n_cols(), n_local_dofs);
4073 }
4074 else
4075 {
4076 std::vector<size_type> &actual_dof_indices = scratch_data->columns;
4077 actual_dof_indices.resize(n_local_dofs);
4078 make_sorted_row_list(local_dof_indices, actual_dof_indices);
4079 const size_type n_actual_dofs = actual_dof_indices.size();
4080 (void)n_actual_dofs;
4081
4082 // additional construct that also takes care of block indices.
4083 std::vector<size_type> &block_starts = scratch_data->block_starts;
4084 block_starts.resize(num_blocks + 1);
4085 internal::AffineConstraints::make_block_starts(sparsity_pattern,
4086 actual_dof_indices,
4087 block_starts);
4088
4089 for (size_type block = 0; block < num_blocks; ++block)
4090 {
4091 const size_type next_block = block_starts[block + 1];
4092 for (size_type i = block_starts[block]; i < next_block; ++i)
4093 {
4094 Assert(i < n_actual_dofs, ExcInternalError());
4095 const size_type row = actual_dof_indices[i];
4096 Assert(row < sparsity_pattern.block(block, 0).n_rows(),
4097 ExcInternalError());
4098 std::vector<size_type>::iterator index_it =
4099 actual_dof_indices.begin();
4100 for (size_type block_col = 0; block_col < num_blocks; ++block_col)
4101 {
4102 const size_type next_block_col = block_starts[block_col + 1];
4103 sparsity_pattern.block(block, block_col)
4104 .add_entries(row,
4105 index_it,
4106 actual_dof_indices.begin() + next_block_col,
4107 true);
4108 index_it = actual_dof_indices.begin() + next_block_col;
4109 }
4110 }
4111 }
4112
4113 for (size_type i = 0; i < n_local_dofs; i++)
4114 if (is_constrained(local_dof_indices[i]))
4115 {
4116 if (keep_constrained_entries == true)
4117 for (size_type j = 0; j < n_local_dofs; j++)
4118 {
4119 sparsity_pattern.add(local_dof_indices[i],
4120 local_dof_indices[j]);
4121 sparsity_pattern.add(local_dof_indices[j],
4122 local_dof_indices[i]);
4123 }
4124 else
4125 sparsity_pattern.add(local_dof_indices[i], local_dof_indices[i]);
4126 }
4127
4128 return;
4129 }
4130
4131 // difficult case with dof_mask, similar to the distribute_local_to_global
4132 // function for block matrices
4133 internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows =
4134 scratch_data->global_rows;
4135 global_rows.reinit(n_local_dofs);
4136 make_sorted_row_list(local_dof_indices, global_rows);
4137 const size_type n_actual_dofs = global_rows.size();
4138
4139 // additional construct that also takes care of block indices.
4140 std::vector<size_type> &block_starts = scratch_data->block_starts;
4141 block_starts.resize(num_blocks + 1);
4142 internal::AffineConstraints::make_block_starts(sparsity_pattern,
4143 global_rows,
4144 block_starts);
4145
4146 std::vector<size_type> &cols = scratch_data->columns;
4147 cols.resize(n_actual_dofs);
4148
4149 // the basic difference to the non-block variant from now onwards is that we
4150 // go through the blocks of the matrix separately.
4151 for (size_type block = 0; block < num_blocks; ++block)
4152 {
4153 const size_type next_block = block_starts[block + 1];
4154 for (size_type i = block_starts[block]; i < next_block; ++i)
4155 {
4156 const size_type row = global_rows.global_row(i);
4157 for (size_type block_col = 0; block_col < num_blocks; ++block_col)
4158 {
4159 const size_type begin_block = block_starts[block_col],
4160 end_block = block_starts[block_col + 1];
4161 std::vector<size_type>::iterator col_ptr = cols.begin();
4162 internal::AffineConstraints::resolve_matrix_row(
4163 global_rows, i, begin_block, end_block, dof_mask, col_ptr);
4164
4165 sparsity_pattern.block(block, block_col)
4166 .add_entries(row, cols.begin(), col_ptr, true);
4167 }
4168 }
4169 }
4170
4171 internal::AffineConstraints::set_sparsity_diagonals(global_rows,
4172 local_dof_indices,
4173 dof_mask,
4174 keep_constrained_entries,
4175 sparsity_pattern);
4176 }
4177
4178 DEAL_II_NAMESPACE_CLOSE
4179
4180 #endif
4181