deal.II version GIT relicensing-2238-gc05b561aad 2024-12-10 20:50:00+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<typename Number::value_type>> 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 -
448 typename Number::value_type(1.)) <
449 100 * std::numeric_limits<double>::epsilon())
450 {
451 current_dof = entries[0].first;
452 goto no_constraint;
453 }
454
455 constraint_indicator.push_back(constraint_iterator);
456 constraint_indicator.back().second =
457 constraint_values.insert_entries(entries);
458
459 // reset constraint iterator for next round
460 constraint_iterator.first = 0;
461
462 if (n_entries > 0)
463 {
464 const std::vector<types::global_dof_index>
465 &constraint_indices = constraint_values.constraint_indices;
466 for (unsigned int j = 0; j < n_entries; ++j)
467 {
468 dof_indices.push_back(
469 global_to_local(constraint_indices[j]));
470 }
471 }
472 }
473 else
474 {
475 no_constraint:
476 dof_indices.push_back(global_to_local(current_dof));
477
478 // make sure constraint_iterator.first is always within the
479 // bounds of unsigned short
480 Assert(constraint_iterator.first <
481 (1 << (8 * sizeof(unsigned short))) - 1,
483 constraint_iterator.first++;
484 }
485 }
486 }
487
488
489
490 template <int dim, typename Number, typename IndexType>
491 inline void
493 const unsigned int cell_no,
494 const std::vector<types::global_dof_index> &local_dof_indices_lex,
495 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
496 {
497 const auto global_to_local =
498 [&](const types::global_dof_index global_index) -> IndexType {
499 if (partitioner)
500 return partitioner->global_to_local(global_index);
501 else
502 {
503 if (local_range.first <= global_index &&
504 global_index < local_range.second)
505 return global_index - local_range.first;
506 else
507 return global_index + (local_range.second - local_range.first);
508 }
509 };
510
511 std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
512
513 auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no];
514 auto &dof_indices = this->dof_indices_per_cell[cell_no];
515
516 for (const auto current_dof : local_dof_indices_lex)
517 {
518 // dof is constrained
519 if (current_dof == numbers::invalid_dof_index)
520 {
521 const std::vector<
522 std::pair<types::global_dof_index, typename Number::value_type>>
523 entries = {};
524
525 constraint_indicator.push_back(constraint_iterator);
526 constraint_indicator.back().second =
527 constraint_values.insert_entries(entries);
528
529 constraint_iterator.first = 0;
530 }
531 else
532 {
533 dof_indices.push_back(global_to_local(current_dof));
534
535 // make sure constraint_iterator.first is always within the
536 // bounds of unsigned short
537 Assert(constraint_iterator.first <
538 (1 << (8 * sizeof(unsigned short))) - 1,
540 constraint_iterator.first++;
541 }
542 }
543 }
544
545
546
547 template <int dim, typename Number, typename IndexType>
548 inline void
550 {
552 std::pair<types::global_dof_index, types::global_dof_index>{0,
553 0}),
555
556 this->dof_indices = {};
557 this->plain_dof_indices = {};
558 this->constraint_indicator = {};
559
560 this->row_starts = {};
561 this->row_starts.emplace_back(0, 0);
562
563 if (this->plain_dof_indices_per_cell.empty() == false)
564 {
565 this->row_starts_plain_indices = {};
566 this->row_starts_plain_indices.emplace_back(0);
567 }
568
569 for (unsigned int i = 0; i < this->dof_indices_per_cell.size(); ++i)
570 {
571 this->dof_indices.insert(this->dof_indices.end(),
572 this->dof_indices_per_cell[i].begin(),
573 this->dof_indices_per_cell[i].end());
574 this->constraint_indicator.insert(
575 this->constraint_indicator.end(),
576 constraint_indicator_per_cell[i].begin(),
577 constraint_indicator_per_cell[i].end());
578
579 this->row_starts.emplace_back(this->dof_indices.size(),
580 this->constraint_indicator.size());
581
582 if (this->plain_dof_indices_per_cell.empty() == false)
583 {
584 this->plain_dof_indices.insert(
585 this->plain_dof_indices.end(),
586 this->plain_dof_indices_per_cell[i].begin(),
587 this->plain_dof_indices_per_cell[i].end());
588
589 this->row_starts_plain_indices.emplace_back(
590 this->plain_dof_indices.size());
591 }
592 }
593
594 std::vector<const std::vector<double> *> constraints(
595 constraint_values.constraints.size());
596 unsigned int length = 0;
597 for (const auto &it : constraint_values.constraints)
598 {
599 AssertIndexRange(it.second, constraints.size());
600 constraints[it.second] = &it.first;
601 length += it.first.size();
602 }
603
604 constraint_pool_data.clear();
605 constraint_pool_data.reserve(length);
606 constraint_pool_row_index.reserve(constraint_values.constraints.size() +
607 1);
608 constraint_pool_row_index.resize(1, 0);
609
610 for (const auto &constraint : constraints)
611 {
612 Assert(constraint != nullptr, ExcInternalError());
613 constraint_pool_data.insert(constraint_pool_data.end(),
614 constraint->begin(),
615 constraint->end());
616 constraint_pool_row_index.push_back(constraint_pool_data.size());
617 }
618
619 AssertDimension(constraint_pool_data.size(), length);
620
621 this->dof_indices_per_cell.clear();
622 constraint_indicator_per_cell.clear();
623
624 if (hanging_nodes &&
625 std::all_of(hanging_node_constraint_masks.begin(),
626 hanging_node_constraint_masks.end(),
627 [](const auto i) {
628 return i == unconstrained_compressed_constraint_kind;
629 }))
630 hanging_node_constraint_masks.clear();
631 }
632
633
634 template <int dim, typename Number, typename IndexType>
635 inline std::shared_ptr<const Utilities::MPI::Partitioner>
637 {
638 this->dof_indices.clear();
639 this->plain_dof_indices.clear();
640 this->constraint_indicator.clear();
641
642 this->row_starts.clear();
643 this->row_starts.reserve(this->dof_indices_per_cell.size());
644 this->row_starts.emplace_back(0, 0);
645
646 if (this->plain_dof_indices_per_cell.empty() == false)
647 {
648 this->row_starts_plain_indices.clear();
649 this->row_starts_plain_indices.reserve(
650 this->dof_indices_per_cell.size());
651 this->row_starts_plain_indices.emplace_back(0);
652 }
653
654 std::vector<types::global_dof_index> ghost_dofs;
655 std::pair<unsigned int, unsigned int> counts = {0, 0};
656
657 for (unsigned int i = 0; i < this->dof_indices_per_cell.size(); ++i)
658 {
659 counts.first += this->dof_indices_per_cell[i].size();
660
661 for (const auto &j : this->dof_indices_per_cell[i])
662 if (j >= (local_range.second - local_range.first))
663 ghost_dofs.push_back(j -
664 (local_range.second - local_range.first));
665
666 if (this->plain_dof_indices_per_cell.empty() == false)
667 {
668 counts.second += this->plain_dof_indices_per_cell[i].size();
669
670 for (const auto &j : this->plain_dof_indices_per_cell[i])
671 if (j >= (local_range.second - local_range.first))
672 ghost_dofs.push_back(
673 j - (local_range.second - local_range.first));
674 }
675 }
676
677 std::sort(ghost_dofs.begin(), ghost_dofs.end());
678 ghost_dofs.erase(std::unique(ghost_dofs.begin(), ghost_dofs.end()),
679 ghost_dofs.end());
680
681 IndexSet locally_relevant_dofs(locally_owned_indices.size());
682 locally_relevant_dofs.add_indices(ghost_dofs.begin(), ghost_dofs.end());
683
684 const auto partitioner =
685 std::make_shared<Utilities::MPI::Partitioner>(locally_owned_indices,
686 locally_relevant_dofs,
687 comm);
688
689 this->dof_indices.reserve(counts.first);
690 this->plain_dof_indices.reserve(counts.second);
691
692 for (unsigned int i = 0; i < this->dof_indices_per_cell.size(); ++i)
693 {
694 for (const auto &j : this->dof_indices_per_cell[i])
695 if (j < (local_range.second - local_range.first))
696 this->dof_indices.push_back(j);
697 else
698 this->dof_indices.push_back(partitioner->global_to_local(
699 j - (local_range.second - local_range.first)));
700
701 this->constraint_indicator.insert(
702 this->constraint_indicator.end(),
703 constraint_indicator_per_cell[i].begin(),
704 constraint_indicator_per_cell[i].end());
705
706 this->row_starts.emplace_back(this->dof_indices.size(),
707 this->constraint_indicator.size());
708
709 if (this->plain_dof_indices_per_cell.empty() == false)
710 {
711 for (const auto &j : this->plain_dof_indices_per_cell[i])
712 if (j < (local_range.second - local_range.first))
713 this->plain_dof_indices.push_back(j);
714 else
715 this->plain_dof_indices.push_back(
716 partitioner->global_to_local(
717 j - (local_range.second - local_range.first)));
718
719 this->row_starts_plain_indices.emplace_back(
720 this->plain_dof_indices.size());
721 }
722 }
723
724 std::vector<const std::vector<double> *> constraints(
725 constraint_values.constraints.size());
726 unsigned int length = 0;
727 for (const auto &it : constraint_values.constraints)
728 {
729 AssertIndexRange(it.second, constraints.size());
730 constraints[it.second] = &it.first;
731 length += it.first.size();
732 }
733
734 constraint_pool_data.clear();
735 constraint_pool_data.reserve(length);
736 constraint_pool_row_index.reserve(constraint_values.constraints.size() +
737 1);
738 constraint_pool_row_index.resize(1, 0);
739
740 for (const auto &constraint : constraints)
741 {
742 Assert(constraint != nullptr, ExcInternalError());
743 constraint_pool_data.insert(constraint_pool_data.end(),
744 constraint->begin(),
745 constraint->end());
746 constraint_pool_row_index.push_back(constraint_pool_data.size());
747 }
748
749 AssertDimension(constraint_pool_data.size(), length);
750
751 this->dof_indices_per_cell.clear();
752 constraint_indicator_per_cell.clear();
753
754 if (hanging_nodes &&
755 std::all_of(hanging_node_constraint_masks.begin(),
756 hanging_node_constraint_masks.end(),
757 [](const auto i) {
758 return i == unconstrained_compressed_constraint_kind;
759 }))
760 hanging_node_constraint_masks.clear();
761
762 return partitioner;
763 }
764
765
766
767 template <int dim, typename Number, typename IndexType>
768 template <typename T, typename VectorType>
769 inline void
771 const T &operation,
772 VectorType &global_vector,
773 Number *local_vector,
774 const unsigned int first_cell,
775 const unsigned int n_cells,
776 const unsigned int n_dofs_per_cell,
777 const bool apply_constraints) const
778 {
779 if ((row_starts_plain_indices.empty() == false) &&
780 (apply_constraints == false))
781 {
782 for (unsigned int v = 0; v < n_cells; ++v)
783 {
784 const unsigned int cell_index = first_cell + v;
785 const unsigned int *dof_indices =
786 this->plain_dof_indices.data() +
787 this->row_starts_plain_indices[cell_index];
788
789 for (unsigned int i = 0; i < n_dofs_per_cell; ++dof_indices, ++i)
790 operation.process_dof(*dof_indices,
791 global_vector,
792 local_vector[i][v]);
793 }
794
795 return;
796 }
797
798 for (unsigned int v = 0; v < n_cells; ++v)
799 {
800 const unsigned int cell_index = first_cell + v;
801 const unsigned int *dof_indices =
802 this->dof_indices.data() + this->row_starts[cell_index].first;
803 unsigned int index_indicators = this->row_starts[cell_index].second;
804 unsigned int next_index_indicators =
805 this->row_starts[cell_index + 1].second;
806
807 unsigned int ind_local = 0;
808 for (; index_indicators != next_index_indicators; ++index_indicators)
809 {
810 const std::pair<unsigned short, unsigned short> indicator =
811 this->constraint_indicator[index_indicators];
812
813 // run through values up to next constraint
814 for (unsigned int j = 0; j < indicator.first; ++j)
815 operation.process_dof(dof_indices[j],
816 global_vector,
817 local_vector[ind_local + j][v]);
818
819 ind_local += indicator.first;
820 dof_indices += indicator.first;
821
822 // constrained case: build the local value as a linear
823 // combination of the global value according to constraints
824 typename Number::value_type value;
825 operation.pre_constraints(local_vector[ind_local][v], value);
826
827 const typename Number::value_type *data_val =
828 this->constraint_pool_begin(indicator.second);
829 const typename Number::value_type *end_pool =
830 this->constraint_pool_end(indicator.second);
831 for (; data_val != end_pool; ++data_val, ++dof_indices)
832 operation.process_constraint(*dof_indices,
833 *data_val,
834 global_vector,
835 value);
836
837 operation.post_constraints(value, local_vector[ind_local][v]);
838 ++ind_local;
839 }
840
841 AssertIndexRange(ind_local, n_dofs_per_cell + 1);
842
843 for (; ind_local < n_dofs_per_cell; ++dof_indices, ++ind_local)
844 operation.process_dof(*dof_indices,
845 global_vector,
846 local_vector[ind_local][v]);
847 }
848 }
849
850
851
852 template <int dim, typename Number, typename IndexType>
853 inline void
855 const unsigned int first_cell,
856 const unsigned int n_lanes_filled,
857 const bool transpose,
858 AlignedVector<Number> &evaluation_data_coarse) const
859 {
860 if (hanging_node_constraint_masks.empty())
861 return;
862
864 Number::size()>
865 constraint_mask;
866
867 bool hn_available = false;
868
869 for (unsigned int v = 0; v < n_lanes_filled; ++v)
870 {
871 const auto mask = hanging_node_constraint_masks[first_cell + v];
872
873 constraint_mask[v] = mask;
874
876 }
877
878 if (hn_available == true)
879 {
880 std::fill(constraint_mask.begin() + n_lanes_filled,
881 constraint_mask.end(),
883
884 for (unsigned int i = 1; i < n_lanes_filled; ++i)
885 AssertDimension(active_fe_indices[first_cell],
886 active_fe_indices[first_cell + i]);
887
888 const auto &shape_info = shape_infos[active_fe_indices[first_cell]];
889
891 dim,
892 typename Number::value_type,
893 Number>::apply(shape_info.n_components,
894 shape_info.data.front().fe_degree,
895 shape_info,
896 transpose,
897 constraint_mask,
898 evaluation_data_coarse.begin());
899 }
900 }
901
902
903
904 template <int dim, typename Number, typename IndexType>
905 inline const typename Number::value_type *
907 const unsigned int row) const
908 {
909 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
910 return constraint_pool_data.empty() ?
911 nullptr :
912 constraint_pool_data.data() + constraint_pool_row_index[row];
913 }
914
915
916
917 template <int dim, typename Number, typename IndexType>
918 inline const typename Number::value_type *
920 const unsigned int row) const
921 {
922 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
923 return constraint_pool_data.empty() ?
924 nullptr :
925 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
926 }
927
928
929
930 template <int dim, typename Number, typename IndexType>
931 inline std::size_t
933 {
934 std::size_t size = 0;
935
936 size += MemoryConsumption::memory_consumption(constraint_values);
937 size += MemoryConsumption::memory_consumption(dof_indices_per_cell);
938 size += MemoryConsumption::memory_consumption(plain_dof_indices_per_cell);
939 size +=
940 MemoryConsumption::memory_consumption(constraint_indicator_per_cell);
941
942 if (hanging_nodes)
944
945 size += MemoryConsumption::memory_consumption(lexicographic_numbering);
947 size += MemoryConsumption::memory_consumption(constraint_indicator);
949 size += MemoryConsumption::memory_consumption(plain_dof_indices);
950 size += MemoryConsumption::memory_consumption(row_starts_plain_indices);
951 size += MemoryConsumption::memory_consumption(constraint_pool_data);
952 size += MemoryConsumption::memory_consumption(constraint_pool_row_index);
954 size +=
955 MemoryConsumption::memory_consumption(hanging_node_constraint_masks);
956 size += MemoryConsumption::memory_consumption(active_fe_indices);
957
958 return size;
959 }
960
961
962
963 template <typename Number>
964 inline std::size_t
966 {
967 std::size_t size = 0;
968
969 size += MemoryConsumption::memory_consumption(constraint_entries);
970 size += MemoryConsumption::memory_consumption(constraint_indices);
971
972 // TODO: map does not have memory_consumption()
973 // size += MemoryConsumption::memory_consumption(constraints);
974
975 return size;
976 }
977
978 } // namespace MatrixFreeFunctions
979} // namespace internal
980
982
983#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:1926
size_type n_elements() const
Definition index_set.h:1934
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1982
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1831
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< ShapeInfo< typename Number::value_type > > shape_infos
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< std::pair< unsigned int, unsigned int > > row_starts
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition grid_out.cc:4630
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
std::pair< types::global_dof_index, types::global_dof_index > local_range
Definition mpi.cc:815
std::size_t size
Definition mpi.cc:734
const MPI_Comm comm
Definition mpi.cc:913
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
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