Reference documentation for deal.II version GIT relicensing-31-gf283a0fc39 2024-02-29 06:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
constraint_info.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_constraint_info_h
17#define dealii_matrix_free_constraint_info_h
18
19
20#include <deal.II/base/config.h>
21
23
28
29#include <limits>
30
32
33namespace internal
34{
35 namespace MatrixFreeFunctions
36 {
41 template <typename Number>
43 {
45
53 template <typename number2>
54 unsigned short
56 const std::vector<std::pair<types::global_dof_index, number2>>
57 &entries);
58
62 std::size_t
64
65 std::vector<std::pair<types::global_dof_index, double>>
67 std::vector<types::global_dof_index> constraint_indices;
68
69 std::pair<std::vector<Number>, types::global_dof_index> next_constraint;
70 std::map<std::vector<Number>,
74 };
75
76
77
83 template <int dim, typename Number, typename IndexType = unsigned int>
85 {
86 public:
91
96 void
98
103 void
104 reinit(const DoFHandler<dim> &dof_handler,
105 const unsigned int n_cells,
106 const bool use_fast_hanging_node_algorithm = true);
107
108 void
110 const unsigned int cell_no,
111 const unsigned int mg_level,
113 const ::AffineConstraints<typename Number::value_type>
114 &constraints,
115 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
116
120 void
121 reinit(const unsigned int n_cells);
122
123 void
125 const unsigned int cell_no,
126 const std::vector<types::global_dof_index> &dof_indices,
127 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
128
129 void
131
132 std::shared_ptr<const Utilities::MPI::Partitioner>
134
135 template <typename T, typename VectorType>
136 void
137 read_write_operation(const T &operation,
138 VectorType &global_vector,
139 Number *local_vector,
140 const unsigned int first_cell,
141 const unsigned int n_cells,
142 const unsigned int n_dofs_per_cell,
143 const bool apply_constraints) const;
144
145 void
147 const unsigned int first_cell,
148 const unsigned int n_lanes_filled,
149 const bool transpose,
150 AlignedVector<Number> &evaluation_data_coarse) const;
151
155 std::size_t
157
158 private:
159 // for setup
161 std::vector<std::vector<unsigned int>> dof_indices_per_cell;
162 std::vector<std::vector<unsigned int>> plain_dof_indices_per_cell;
163 std::vector<std::vector<std::pair<unsigned short, unsigned short>>>
165
166 std::unique_ptr<HangingNodes<dim>> hanging_nodes;
167 std::vector<std::vector<unsigned int>> lexicographic_numbering;
168
169 std::vector<types::global_dof_index> local_dof_indices;
170 std::vector<types::global_dof_index> local_dof_indices_lex;
171 std::vector<ConstraintKinds> mask;
172
174 std::pair<types::global_dof_index, types::global_dof_index> local_range;
175
176 public:
177 // for read_write_operation()
178 std::vector<unsigned int> dof_indices;
179 std::vector<std::pair<unsigned short, unsigned short>>
181 std::vector<std::pair<unsigned int, unsigned int>> row_starts;
182
183 std::vector<unsigned int> plain_dof_indices;
184 std::vector<unsigned int> row_starts_plain_indices;
185
186 // for constraint_pool_begin/end()
187 std::vector<typename Number::value_type> constraint_pool_data;
188 std::vector<unsigned int> constraint_pool_row_index;
189
190 std::vector<ShapeInfo<Number>> shape_infos;
191 std::vector<compressed_constraint_kind> hanging_node_constraint_masks;
192 std::vector<unsigned int> active_fe_indices;
193
194 private:
195 inline const typename Number::value_type *
196 constraint_pool_begin(const unsigned int row) const;
197
198 inline const typename Number::value_type *
199 constraint_pool_end(const unsigned int row) const;
200 };
201
202
203
204 // ------------------------- inline functions --------------------------
205
206 // NOLINTNEXTLINE(modernize-use-equals-default)
207 template <typename Number>
209 : constraints(FloatingPointComparator<Number>(
210 1. * std::numeric_limits<double>::epsilon() * 1024.))
211 {}
212
213
214
215 template <typename Number>
216 template <typename number2>
217 unsigned short
219 const std::vector<std::pair<types::global_dof_index, number2>> &entries)
220 {
221 next_constraint.first.resize(entries.size());
222 if (entries.size() > 0)
223 {
224 constraint_indices.resize(entries.size());
225 // Use assign so that values for nonmatching Number / number2 are
226 // converted:
227 constraint_entries.assign(entries.begin(), entries.end());
228 std::sort(constraint_entries.begin(),
229 constraint_entries.end(),
230 [](const std::pair<types::global_dof_index, double> &p1,
231 const std::pair<types::global_dof_index, double> &p2) {
232 return p1.second < p2.second;
233 });
234 for (types::global_dof_index j = 0; j < constraint_entries.size();
235 j++)
236 {
237 // copy the indices of the constraint entries after sorting.
238 constraint_indices[j] = constraint_entries[j].first;
239
240 // one_constraint takes the weights of the constraint
241 next_constraint.first[j] = constraint_entries[j].second;
242 }
243 }
244
245 // check whether or not constraint is already in pool. the initial
246 // implementation computed a hash value based on the truncated array (to
247 // given accuracy around 1e-13) in order to easily detect different
248 // arrays and then made a fine-grained check when the hash values were
249 // equal. this was quite lengthy and now we use a std::map with a
250 // user-defined comparator to compare floating point arrays to a
251 // tolerance 1e-13.
253 const auto position = constraints.find(next_constraint.first);
254 if (position != constraints.end())
255 insert_position = position->second;
256 else
257 {
258 next_constraint.second = constraints.size();
259 constraints.insert(next_constraint);
260 insert_position = next_constraint.second;
261 }
262
263 // we want to store the result as a short variable, so we have to make
264 // sure that the result does not exceed the limits when casting.
265 Assert(insert_position < (1 << (8 * sizeof(unsigned short))),
267 return static_cast<unsigned short>(insert_position);
268 }
269
270
271
272 template <int dim, typename Number, typename IndexType>
276
277
278
279 template <int dim, typename Number, typename IndexType>
280 void
282 const IndexSet &locally_owned_indices)
283 {
284 this->locally_owned_indices = locally_owned_indices;
285
286 if (locally_owned_indices.is_empty())
287 local_range = {0, 0};
288 else
289 local_range = {locally_owned_indices.nth_index_in_set(0),
290 locally_owned_indices.nth_index_in_set(0) +
291 locally_owned_indices.n_elements()};
292 }
293
294
295
296 template <int dim, typename Number, typename IndexType>
297 inline void
299 const DoFHandler<dim> &dof_handler,
300 const unsigned int n_cells,
301 const bool use_fast_hanging_node_algorithm)
302 {
303 this->dof_indices_per_cell.resize(n_cells);
304 this->plain_dof_indices_per_cell.resize(n_cells);
305 this->constraint_indicator_per_cell.resize(n_cells);
306
307 // note: has_hanging_nodes() is a global operatrion
308 const bool has_hanging_nodes =
309 dof_handler.get_triangulation().has_hanging_nodes();
310
311 if (use_fast_hanging_node_algorithm && has_hanging_nodes)
312 {
313 hanging_nodes = std::make_unique<HangingNodes<dim>>(
314 dof_handler.get_triangulation());
315
316 hanging_node_constraint_masks.resize(n_cells);
317 }
318
319 const auto &fes = dof_handler.get_fe_collection();
320 lexicographic_numbering.resize(fes.size());
321 shape_infos.resize(fes.size());
322
323 for (unsigned int i = 0; i < fes.size(); ++i)
324 {
325 if (fes[i].reference_cell().is_hyper_cube())
326 {
327 const Quadrature<1> dummy_quadrature(
328 std::vector<Point<1>>(1, Point<1>()));
329 shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
330 }
331 else
332 {
333 const auto dummy_quadrature =
334 fes[i].reference_cell().template get_gauss_type_quadrature<dim>(
335 1);
336 shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
337 }
338
339 lexicographic_numbering[i] = shape_infos[i].lexicographic_numbering;
340 }
341 active_fe_indices.resize(n_cells);
342 }
343
344
345
346 template <int dim, typename Number, typename IndexType>
347 inline void
349 {
350 this->dof_indices_per_cell.resize(n_cells);
351 this->plain_dof_indices_per_cell.resize(0);
352 this->constraint_indicator_per_cell.resize(n_cells);
353 }
354
355
356
357 template <int dim, typename Number, typename IndexType>
358 inline void
360 const unsigned int cell_no,
361 const unsigned int mg_level,
363 const ::AffineConstraints<typename Number::value_type> &constraints,
364 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
365 {
366 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
367 local_dof_indices_lex.resize(cell->get_fe().n_dofs_per_cell());
368
369 if (mg_level == numbers::invalid_unsigned_int)
370 cell->get_dof_indices(local_dof_indices);
371 else
372 cell->get_mg_dof_indices(local_dof_indices);
373
374 {
375 AssertIndexRange(cell->active_fe_index(), shape_infos.size());
376
377 const auto &lexicographic_numbering =
378 shape_infos[cell->active_fe_index()].lexicographic_numbering;
379
380 AssertDimension(lexicographic_numbering.size(),
381 local_dof_indices.size());
382
383 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
384 local_dof_indices_lex[i] =
385 local_dof_indices[lexicographic_numbering[i]];
386 }
387
388 std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
389
390 AssertIndexRange(cell_no, this->constraint_indicator_per_cell.size());
391 AssertIndexRange(cell_no, this->dof_indices_per_cell.size());
392 AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size());
393
394 auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no];
395 auto &dof_indices = this->dof_indices_per_cell[cell_no];
396 auto &plain_dof_indices = this->plain_dof_indices_per_cell[cell_no];
397
398 AssertDimension(constraint_indicator_per_cell[cell_no].size(), 0);
399 AssertDimension(this->dof_indices_per_cell[cell_no].size(), 0);
400 AssertDimension(this->plain_dof_indices_per_cell[cell_no].size(), 0);
401
402 const auto global_to_local =
403 [&](const types::global_dof_index global_index) -> IndexType {
404 if (partitioner)
405 return partitioner->global_to_local(global_index);
406 else
407 {
408 if (local_range.first <= global_index &&
409 global_index < local_range.second)
410 return global_index - local_range.first;
411 else
412 return global_index + (local_range.second - local_range.first);
413 }
414 };
415
416 // plain indices
417 plain_dof_indices.resize(local_dof_indices_lex.size());
418 for (unsigned int i = 0; i < local_dof_indices_lex.size(); ++i)
419 plain_dof_indices[i] = global_to_local(local_dof_indices_lex[i]);
420
421 if (hanging_nodes)
422 {
423 AssertIndexRange(cell_no, this->hanging_node_constraint_masks.size());
424 AssertIndexRange(cell_no, this->active_fe_indices.size());
425
426 mask.assign(
427 cell->get_fe().n_components(),
429 hanging_nodes->setup_constraints(
430 cell, {}, lexicographic_numbering, local_dof_indices_lex, mask);
431
432 hanging_node_constraint_masks[cell_no] = compress(mask[0], dim);
433 active_fe_indices[cell_no] = cell->active_fe_index();
434 }
435
436 for (auto current_dof : local_dof_indices_lex)
437 {
438 const auto *entries_ptr =
439 constraints.get_constraint_entries(current_dof);
440
441 // dof is constrained
442 if (entries_ptr != nullptr)
443 {
444 const auto &entries = *entries_ptr;
445 const types::global_dof_index n_entries = entries.size();
446 if (n_entries == 1 &&
447 std::abs(entries[0].second - 1.) <
448 100 * std::numeric_limits<double>::epsilon())
449 {
450 current_dof = entries[0].first;
451 goto no_constraint;
452 }
453
454 constraint_indicator.push_back(constraint_iterator);
455 constraint_indicator.back().second =
456 constraint_values.insert_entries(entries);
457
458 // reset constraint iterator for next round
459 constraint_iterator.first = 0;
460
461 if (n_entries > 0)
462 {
463 const std::vector<types::global_dof_index>
464 &constraint_indices = constraint_values.constraint_indices;
465 for (unsigned int j = 0; j < n_entries; ++j)
466 {
467 dof_indices.push_back(
468 global_to_local(constraint_indices[j]));
469 }
470 }
471 }
472 else
473 {
474 no_constraint:
475 dof_indices.push_back(global_to_local(current_dof));
476
477 // make sure constraint_iterator.first is always within the
478 // bounds of unsigned short
479 Assert(constraint_iterator.first <
480 (1 << (8 * sizeof(unsigned short))) - 1,
482 constraint_iterator.first++;
483 }
484 }
485 }
486
487
488
489 template <int dim, typename Number, typename IndexType>
490 inline void
492 const unsigned int cell_no,
493 const std::vector<types::global_dof_index> &local_dof_indices_lex,
494 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
495 {
496 const auto global_to_local =
497 [&](const types::global_dof_index global_index) -> IndexType {
498 if (partitioner)
499 return partitioner->global_to_local(global_index);
500 else
501 {
502 if (local_range.first <= global_index &&
503 global_index < local_range.second)
504 return global_index - local_range.first;
505 else
506 return global_index + (local_range.second - local_range.first);
507 }
508 };
509
510 std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
511
512 auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no];
513 auto &dof_indices = this->dof_indices_per_cell[cell_no];
514
515 for (const auto current_dof : local_dof_indices_lex)
516 {
517 // dof is constrained
518 if (current_dof == numbers::invalid_dof_index)
519 {
520 const std::vector<
521 std::pair<types::global_dof_index, typename Number::value_type>>
522 entries = {};
523
524 constraint_indicator.push_back(constraint_iterator);
525 constraint_indicator.back().second =
526 constraint_values.insert_entries(entries);
527
528 constraint_iterator.first = 0;
529 }
530 else
531 {
532 dof_indices.push_back(global_to_local(current_dof));
533
534 // make sure constraint_iterator.first is always within the
535 // bounds of unsigned short
536 Assert(constraint_iterator.first <
537 (1 << (8 * sizeof(unsigned short))) - 1,
539 constraint_iterator.first++;
540 }
541 }
542 }
543
544
545
546 template <int dim, typename Number, typename IndexType>
547 inline void
549 {
550 Assert((local_range ==
551 std::pair<types::global_dof_index, types::global_dof_index>{0,
552 0}),
554
555 this->dof_indices = {};
556 this->plain_dof_indices = {};
557 this->constraint_indicator = {};
558
559 this->row_starts = {};
560 this->row_starts.emplace_back(0, 0);
561
562 if (this->plain_dof_indices_per_cell.empty() == false)
563 {
564 this->row_starts_plain_indices = {};
565 this->row_starts_plain_indices.emplace_back(0);
566 }
567
568 for (unsigned int i = 0; i < this->dof_indices_per_cell.size(); ++i)
569 {
570 this->dof_indices.insert(this->dof_indices.end(),
571 this->dof_indices_per_cell[i].begin(),
572 this->dof_indices_per_cell[i].end());
573 this->constraint_indicator.insert(
574 this->constraint_indicator.end(),
575 constraint_indicator_per_cell[i].begin(),
576 constraint_indicator_per_cell[i].end());
577
578 this->row_starts.emplace_back(this->dof_indices.size(),
579 this->constraint_indicator.size());
580
581 if (this->plain_dof_indices_per_cell.empty() == false)
582 {
583 this->plain_dof_indices.insert(
584 this->plain_dof_indices.end(),
585 this->plain_dof_indices_per_cell[i].begin(),
586 this->plain_dof_indices_per_cell[i].end());
587
588 this->row_starts_plain_indices.emplace_back(
589 this->plain_dof_indices.size());
590 }
591 }
592
593 std::vector<const std::vector<double> *> constraints(
594 constraint_values.constraints.size());
595 unsigned int length = 0;
596 for (const auto &it : constraint_values.constraints)
597 {
598 AssertIndexRange(it.second, constraints.size());
599 constraints[it.second] = &it.first;
600 length += it.first.size();
601 }
602
603 constraint_pool_data.clear();
604 constraint_pool_data.reserve(length);
605 constraint_pool_row_index.reserve(constraint_values.constraints.size() +
606 1);
607 constraint_pool_row_index.resize(1, 0);
608
609 for (const auto &constraint : constraints)
610 {
611 Assert(constraint != nullptr, ExcInternalError());
612 constraint_pool_data.insert(constraint_pool_data.end(),
613 constraint->begin(),
614 constraint->end());
615 constraint_pool_row_index.push_back(constraint_pool_data.size());
616 }
617
618 AssertDimension(constraint_pool_data.size(), length);
619
620 this->dof_indices_per_cell.clear();
621 constraint_indicator_per_cell.clear();
622
623 if (hanging_nodes &&
624 std::all_of(hanging_node_constraint_masks.begin(),
625 hanging_node_constraint_masks.end(),
626 [](const auto i) {
627 return i == unconstrained_compressed_constraint_kind;
628 }))
629 hanging_node_constraint_masks.clear();
630 }
631
632
633 template <int dim, typename Number, typename IndexType>
634 inline std::shared_ptr<const Utilities::MPI::Partitioner>
636 {
637 this->dof_indices.clear();
638 this->plain_dof_indices.clear();
639 this->constraint_indicator.clear();
640
641 this->row_starts.clear();
642 this->row_starts.reserve(this->dof_indices_per_cell.size());
643 this->row_starts.emplace_back(0, 0);
644
645 if (this->plain_dof_indices_per_cell.empty() == false)
646 {
647 this->row_starts_plain_indices.clear();
648 this->row_starts_plain_indices.reserve(
649 this->dof_indices_per_cell.size());
650 this->row_starts_plain_indices.emplace_back(0);
651 }
652
653 std::vector<types::global_dof_index> ghost_dofs;
654 std::pair<unsigned int, unsigned int> counts = {0, 0};
655
656 for (unsigned int i = 0; i < this->dof_indices_per_cell.size(); ++i)
657 {
658 counts.first += this->dof_indices_per_cell[i].size();
659
660 for (const auto &j : this->dof_indices_per_cell[i])
661 if (j >= (local_range.second - local_range.first))
662 ghost_dofs.push_back(j -
663 (local_range.second - local_range.first));
664
665 if (this->plain_dof_indices_per_cell.empty() == false)
666 {
667 counts.second += this->plain_dof_indices_per_cell[i].size();
668
669 for (const auto &j : this->plain_dof_indices_per_cell[i])
670 if (j >= (local_range.second - local_range.first))
671 ghost_dofs.push_back(
672 j - (local_range.second - local_range.first));
673 }
674 }
675
676 std::sort(ghost_dofs.begin(), ghost_dofs.end());
677 ghost_dofs.erase(std::unique(ghost_dofs.begin(), ghost_dofs.end()),
678 ghost_dofs.end());
679
680 IndexSet locally_relevant_dofs(locally_owned_indices.size());
681 locally_relevant_dofs.add_indices(ghost_dofs.begin(), ghost_dofs.end());
682
683 const auto partitioner =
684 std::make_shared<Utilities::MPI::Partitioner>(locally_owned_indices,
685 locally_relevant_dofs,
686 comm);
687
688 this->dof_indices.reserve(counts.first);
689 this->plain_dof_indices.reserve(counts.second);
690
691 for (unsigned int i = 0; i < this->dof_indices_per_cell.size(); ++i)
692 {
693 for (const auto &j : this->dof_indices_per_cell[i])
694 if (j < (local_range.second - local_range.first))
695 this->dof_indices.push_back(j);
696 else
697 this->dof_indices.push_back(partitioner->global_to_local(
698 j - (local_range.second - local_range.first)));
699
700 this->constraint_indicator.insert(
701 this->constraint_indicator.end(),
702 constraint_indicator_per_cell[i].begin(),
703 constraint_indicator_per_cell[i].end());
704
705 this->row_starts.emplace_back(this->dof_indices.size(),
706 this->constraint_indicator.size());
707
708 if (this->plain_dof_indices_per_cell.empty() == false)
709 {
710 for (const auto &j : this->plain_dof_indices_per_cell[i])
711 if (j < (local_range.second - local_range.first))
712 this->plain_dof_indices.push_back(j);
713 else
714 this->plain_dof_indices.push_back(
715 partitioner->global_to_local(
716 j - (local_range.second - local_range.first)));
717
718 this->row_starts_plain_indices.emplace_back(
719 this->plain_dof_indices.size());
720 }
721 }
722
723 std::vector<const std::vector<double> *> constraints(
724 constraint_values.constraints.size());
725 unsigned int length = 0;
726 for (const auto &it : constraint_values.constraints)
727 {
728 AssertIndexRange(it.second, constraints.size());
729 constraints[it.second] = &it.first;
730 length += it.first.size();
731 }
732
733 constraint_pool_data.clear();
734 constraint_pool_data.reserve(length);
735 constraint_pool_row_index.reserve(constraint_values.constraints.size() +
736 1);
737 constraint_pool_row_index.resize(1, 0);
738
739 for (const auto &constraint : constraints)
740 {
741 Assert(constraint != nullptr, ExcInternalError());
742 constraint_pool_data.insert(constraint_pool_data.end(),
743 constraint->begin(),
744 constraint->end());
745 constraint_pool_row_index.push_back(constraint_pool_data.size());
746 }
747
748 AssertDimension(constraint_pool_data.size(), length);
749
750 this->dof_indices_per_cell.clear();
751 constraint_indicator_per_cell.clear();
752
753 if (hanging_nodes &&
754 std::all_of(hanging_node_constraint_masks.begin(),
755 hanging_node_constraint_masks.end(),
756 [](const auto i) {
757 return i == unconstrained_compressed_constraint_kind;
758 }))
759 hanging_node_constraint_masks.clear();
760
761 return partitioner;
762 }
763
764
765
766 template <int dim, typename Number, typename IndexType>
767 template <typename T, typename VectorType>
768 inline void
770 const T &operation,
771 VectorType &global_vector,
772 Number *local_vector,
773 const unsigned int first_cell,
774 const unsigned int n_cells,
775 const unsigned int n_dofs_per_cell,
776 const bool apply_constraints) const
777 {
778 if ((row_starts_plain_indices.empty() == false) &&
779 (apply_constraints == false))
780 {
781 for (unsigned int v = 0; v < n_cells; ++v)
782 {
783 const unsigned int cell_index = first_cell + v;
784 const unsigned int *dof_indices =
785 this->plain_dof_indices.data() +
786 this->row_starts_plain_indices[cell_index];
787
788 for (unsigned int i = 0; i < n_dofs_per_cell; ++dof_indices, ++i)
789 operation.process_dof(*dof_indices,
790 global_vector,
791 local_vector[i][v]);
792 }
793
794 return;
795 }
796
797 for (unsigned int v = 0; v < n_cells; ++v)
798 {
799 const unsigned int cell_index = first_cell + v;
800 const unsigned int *dof_indices =
801 this->dof_indices.data() + this->row_starts[cell_index].first;
802 unsigned int index_indicators = this->row_starts[cell_index].second;
803 unsigned int next_index_indicators =
804 this->row_starts[cell_index + 1].second;
805
806 unsigned int ind_local = 0;
807 for (; index_indicators != next_index_indicators; ++index_indicators)
808 {
809 const std::pair<unsigned short, unsigned short> indicator =
810 this->constraint_indicator[index_indicators];
811
812 // run through values up to next constraint
813 for (unsigned int j = 0; j < indicator.first; ++j)
814 operation.process_dof(dof_indices[j],
815 global_vector,
816 local_vector[ind_local + j][v]);
817
818 ind_local += indicator.first;
819 dof_indices += indicator.first;
820
821 // constrained case: build the local value as a linear
822 // combination of the global value according to constraints
823 typename Number::value_type value;
824 operation.pre_constraints(local_vector[ind_local][v], value);
825
826 const typename Number::value_type *data_val =
827 this->constraint_pool_begin(indicator.second);
828 const typename Number::value_type *end_pool =
829 this->constraint_pool_end(indicator.second);
830 for (; data_val != end_pool; ++data_val, ++dof_indices)
831 operation.process_constraint(*dof_indices,
832 *data_val,
833 global_vector,
834 value);
835
836 operation.post_constraints(value, local_vector[ind_local][v]);
837 ++ind_local;
838 }
839
840 AssertIndexRange(ind_local, n_dofs_per_cell + 1);
841
842 for (; ind_local < n_dofs_per_cell; ++dof_indices, ++ind_local)
843 operation.process_dof(*dof_indices,
844 global_vector,
845 local_vector[ind_local][v]);
846 }
847 }
848
849
850
851 template <int dim, typename Number, typename IndexType>
852 inline void
854 const unsigned int first_cell,
855 const unsigned int n_lanes_filled,
856 const bool transpose,
857 AlignedVector<Number> &evaluation_data_coarse) const
858 {
859 if (hanging_node_constraint_masks.empty())
860 return;
861
863 Number::size()>
864 constraint_mask;
865
866 bool hn_available = false;
867
868 for (unsigned int v = 0; v < n_lanes_filled; ++v)
869 {
870 const auto mask = hanging_node_constraint_masks[first_cell + v];
871
872 constraint_mask[v] = mask;
873
875 }
876
877 if (hn_available == true)
878 {
879 std::fill(constraint_mask.begin() + n_lanes_filled,
880 constraint_mask.end(),
882
883 for (unsigned int i = 1; i < n_lanes_filled; ++i)
884 AssertDimension(active_fe_indices[first_cell],
885 active_fe_indices[first_cell + i]);
886
887 const auto &shape_info = shape_infos[active_fe_indices[first_cell]];
888
890 dim,
891 typename Number::value_type,
892 Number>::apply(shape_info.n_components,
893 shape_info.data.front().fe_degree,
894 shape_info,
895 transpose,
896 constraint_mask,
897 evaluation_data_coarse.begin());
898 }
899 }
900
901
902
903 template <int dim, typename Number, typename IndexType>
904 inline const typename Number::value_type *
906 const unsigned int row) const
907 {
908 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
909 return constraint_pool_data.empty() ?
910 nullptr :
911 constraint_pool_data.data() + constraint_pool_row_index[row];
912 }
913
914
915
916 template <int dim, typename Number, typename IndexType>
917 inline const typename Number::value_type *
919 const unsigned int row) const
920 {
921 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
922 return constraint_pool_data.empty() ?
923 nullptr :
924 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
925 }
926
927
928
929 template <int dim, typename Number, typename IndexType>
930 inline std::size_t
932 {
933 std::size_t size = 0;
934
935 size += MemoryConsumption::memory_consumption(constraint_values);
936 size += MemoryConsumption::memory_consumption(dof_indices_per_cell);
937 size += MemoryConsumption::memory_consumption(plain_dof_indices_per_cell);
938 size +=
939 MemoryConsumption::memory_consumption(constraint_indicator_per_cell);
940
941 if (hanging_nodes)
942 size += MemoryConsumption::memory_consumption(*hanging_nodes);
943
944 size += MemoryConsumption::memory_consumption(lexicographic_numbering);
945 size += MemoryConsumption::memory_consumption(dof_indices);
946 size += MemoryConsumption::memory_consumption(constraint_indicator);
947 size += MemoryConsumption::memory_consumption(row_starts);
948 size += MemoryConsumption::memory_consumption(plain_dof_indices);
949 size += MemoryConsumption::memory_consumption(row_starts_plain_indices);
950 size += MemoryConsumption::memory_consumption(constraint_pool_data);
951 size += MemoryConsumption::memory_consumption(constraint_pool_row_index);
952 size += MemoryConsumption::memory_consumption(shape_infos);
953 size +=
954 MemoryConsumption::memory_consumption(hanging_node_constraint_masks);
955 size += MemoryConsumption::memory_consumption(active_fe_indices);
956
957 return size;
958 }
959
960
961
962 template <typename Number>
963 inline std::size_t
965 {
966 std::size_t size = 0;
967
968 size += MemoryConsumption::memory_consumption(constraint_entries);
969 size += MemoryConsumption::memory_consumption(constraint_indices);
970
971 // TODO: map does not have memory_consumption()
972 // size += MemoryConsumption::memory_consumption(constraints);
973
974 return size;
975 }
976
977 } // namespace MatrixFreeFunctions
978} // namespace internal
979
981
982#endif
iterator begin()
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
bool is_empty() const
Definition index_set.h:1915
size_type n_elements() const
Definition index_set.h:1923
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1971
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1820
Definition point.h:111
virtual bool has_hanging_nodes() const
void apply_hanging_node_constraints(const unsigned int first_cell, const unsigned int n_lanes_filled, const bool transpose, AlignedVector< Number > &evaluation_data_coarse) const
std::vector< std::vector< unsigned int > > lexicographic_numbering
void read_write_operation(const T &operation, VectorType &global_vector, Number *local_vector, const unsigned int first_cell, const unsigned int n_cells, const unsigned int n_dofs_per_cell, const bool apply_constraints) const
std::vector< typename Number::value_type > constraint_pool_data
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
std::vector< std::vector< unsigned int > > dof_indices_per_cell
std::vector< std::vector< unsigned int > > plain_dof_indices_per_cell
const Number::value_type * constraint_pool_begin(const unsigned int row) const
std::vector< types::global_dof_index > local_dof_indices_lex
void set_locally_owned_indices(const IndexSet &locally_owned_indices)
void read_dof_indices(const unsigned int cell_no, const std::vector< types::global_dof_index > &dof_indices, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner)
std::pair< types::global_dof_index, types::global_dof_index > local_range
const Number::value_type * constraint_pool_end(const unsigned int row) const
std::vector< types::global_dof_index > local_dof_indices
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
std::shared_ptr< const Utilities::MPI::Partitioner > finalize(const MPI_Comm comm)
std::vector< unsigned int > constraint_pool_row_index
std::unique_ptr< HangingNodes< dim > > hanging_nodes
void read_dof_indices(const unsigned int cell_no, const unsigned int mg_level, const TriaIterator< DoFCellAccessor< dim, dim, false > > &cell, const ::AffineConstraints< typename Number::value_type > &constraints, const std::shared_ptr< const Utilities::MPI::Partitioner > &partitioner)
std::vector< std::vector< std::pair< unsigned short, unsigned short > > > constraint_indicator_per_cell
void reinit(const unsigned int n_cells)
void reinit(const DoFHandler< dim > &dof_handler, const unsigned int n_cells, const bool use_fast_hanging_node_algorithm=true)
std::vector< unsigned int > row_starts_plain_indices
std::vector< ShapeInfo< Number > > shape_infos
std::vector< std::pair< unsigned int, unsigned int > > row_starts
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition grid_out.cc:4614
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
std::uint8_t compressed_constraint_kind
Definition dof_info.h:86
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:81
*braid_SplitCommworld & comm
std::vector< types::global_dof_index > constraint_indices
std::map< std::vector< Number >, types::global_dof_index, FloatingPointComparator< Number > > constraints
unsigned short insert_entries(const std::vector< std::pair< types::global_dof_index, number2 > > &entries)
std::vector< std::pair< types::global_dof_index, double > > constraint_entries
std::pair< std::vector< Number >, types::global_dof_index > next_constraint