16#ifndef dealii_matrix_free_constraint_info_h
17#define dealii_matrix_free_constraint_info_h
41 template <
typename Number>
53 template <
typename number2>
56 const std::vector<std::pair<types::global_dof_index, number2>>
65 std::vector<std::pair<types::global_dof_index, double>>
70 std::map<std::vector<Number>,
83 template <
int dim,
typename Number,
typename IndexType =
unsigned int>
105 const unsigned int n_cells,
106 const bool use_fast_hanging_node_algorithm =
true);
110 const unsigned int cell_no,
111 const unsigned int mg_level,
113 const ::AffineConstraints<typename Number::value_type>
115 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
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);
132 std::shared_ptr<const Utilities::MPI::Partitioner>
135 template <
typename T,
typename VectorType>
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;
147 const unsigned int first_cell,
148 const unsigned int n_lanes_filled,
163 std::vector<std::vector<std::pair<unsigned short, unsigned short>>>
171 std::vector<ConstraintKinds>
mask;
174 std::pair<types::global_dof_index, types::global_dof_index>
local_range;
179 std::vector<std::pair<unsigned short, unsigned short>>
181 std::vector<std::pair<unsigned int, unsigned int>>
row_starts;
195 inline const typename Number::value_type *
198 inline const typename Number::value_type *
207 template <
typename Number>
210 1. *
std::numeric_limits<double>::epsilon() * 1024.))
215 template <
typename Number>
216 template <
typename number2>
219 const std::vector<std::pair<types::global_dof_index, number2>> &entries)
222 if (entries.size() > 0)
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;
255 insert_position = position->second;
265 Assert(insert_position < (1 << (8 *
sizeof(
unsigned short))),
267 return static_cast<unsigned short>(insert_position);
272 template <
int dim,
typename Number,
typename IndexType>
279 template <
int dim,
typename Number,
typename IndexType>
296 template <
int dim,
typename Number,
typename IndexType>
300 const unsigned int n_cells,
301 const bool use_fast_hanging_node_algorithm)
308 const bool has_hanging_nodes =
311 if (use_fast_hanging_node_algorithm && has_hanging_nodes)
323 for (
unsigned int i = 0; i < fes.size(); ++i)
325 if (fes[i].reference_cell().is_hyper_cube())
329 shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
333 const auto dummy_quadrature =
334 fes[i].reference_cell().template get_gauss_type_quadrature<dim>(
336 shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
346 template <
int dim,
typename Number,
typename IndexType>
357 template <
int dim,
typename Number,
typename IndexType>
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)
378 shape_infos[cell->active_fe_index()].lexicographic_numbering;
383 for (
unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
388 std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
402 const auto global_to_local =
405 return partitioner->global_to_local(global_index);
427 cell->get_fe().n_components(),
438 const auto *entries_ptr =
439 constraints.get_constraint_entries(current_dof);
442 if (entries_ptr !=
nullptr)
444 const auto &entries = *entries_ptr;
446 if (n_entries == 1 &&
448 typename Number::value_type(1.)) <
449 100 * std::numeric_limits<double>::epsilon())
451 current_dof = entries[0].first;
460 constraint_iterator.first = 0;
464 const std::vector<types::global_dof_index>
466 for (
unsigned int j = 0; j < n_entries; ++j)
469 global_to_local(constraint_indices[j]));
476 dof_indices.push_back(global_to_local(current_dof));
480 Assert(constraint_iterator.first <
481 (1 << (8 *
sizeof(
unsigned short))) - 1,
483 constraint_iterator.first++;
490 template <
int dim,
typename Number,
typename IndexType>
493 const unsigned int cell_no,
495 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
497 const auto global_to_local =
500 return partitioner->global_to_local(global_index);
511 std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
522 std::pair<types::global_dof_index, typename Number::value_type>>
529 constraint_iterator.first = 0;
533 dof_indices.push_back(global_to_local(current_dof));
537 Assert(constraint_iterator.first <
538 (1 << (8 *
sizeof(
unsigned short))) - 1,
540 constraint_iterator.first++;
547 template <
int dim,
typename Number,
typename IndexType>
552 std::pair<types::global_dof_index, types::global_dof_index>{0,
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());
579 this->
row_starts.emplace_back(this->dof_indices.size(),
580 this->constraint_indicator.size());
586 this->plain_dof_indices_per_cell[i].begin(),
587 this->plain_dof_indices_per_cell[i].end());
594 std::vector<const std::vector<double> *> constraints(
596 unsigned int length = 0;
600 constraints[it.second] = &it.first;
601 length += it.first.size();
610 for (
const auto &constraint : constraints)
628 return i == unconstrained_compressed_constraint_kind;
634 template <
int dim,
typename Number,
typename IndexType>
635 inline std::shared_ptr<const Utilities::MPI::Partitioner>
654 std::vector<types::global_dof_index> ghost_dofs;
655 std::pair<unsigned int, unsigned int> counts = {0, 0};
663 ghost_dofs.push_back(j -
672 ghost_dofs.push_back(
677 std::sort(ghost_dofs.begin(), ghost_dofs.end());
678 ghost_dofs.erase(std::unique(ghost_dofs.begin(), ghost_dofs.end()),
682 locally_relevant_dofs.
add_indices(ghost_dofs.begin(), ghost_dofs.end());
684 const auto partitioner =
686 locally_relevant_dofs,
698 this->
dof_indices.push_back(partitioner->global_to_local(
707 this->constraint_indicator.size());
716 partitioner->global_to_local(
724 std::vector<const std::vector<double> *> constraints(
726 unsigned int length = 0;
730 constraints[it.second] = &it.first;
731 length += it.first.size();
740 for (
const auto &constraint : constraints)
758 return i == unconstrained_compressed_constraint_kind;
767 template <
int dim,
typename Number,
typename IndexType>
768 template <
typename T,
typename VectorType>
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
780 (apply_constraints ==
false))
782 for (
unsigned int v = 0; v < n_cells; ++v)
784 const unsigned int cell_index = first_cell + v;
789 for (
unsigned int i = 0; i < n_dofs_per_cell; ++
dof_indices, ++i)
798 for (
unsigned int v = 0; v < n_cells; ++v)
800 const unsigned int cell_index = first_cell + v;
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 =
807 unsigned int ind_local = 0;
808 for (; index_indicators != next_index_indicators; ++index_indicators)
810 const std::pair<unsigned short, unsigned short> indicator =
814 for (
unsigned int j = 0; j < indicator.first; ++j)
817 local_vector[ind_local + j][v]);
819 ind_local += indicator.first;
824 typename Number::value_type
value;
825 operation.pre_constraints(local_vector[ind_local][v],
value);
827 const typename Number::value_type *data_val =
829 const typename Number::value_type *end_pool =
831 for (; data_val != end_pool; ++data_val, ++
dof_indices)
837 operation.post_constraints(
value, local_vector[ind_local][v]);
843 for (; ind_local < n_dofs_per_cell; ++
dof_indices, ++ind_local)
846 local_vector[ind_local][v]);
852 template <
int dim,
typename Number,
typename IndexType>
855 const unsigned int first_cell,
856 const unsigned int n_lanes_filled,
867 bool hn_available =
false;
869 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
873 constraint_mask[v] =
mask;
878 if (hn_available ==
true)
880 std::fill(constraint_mask.begin() + n_lanes_filled,
881 constraint_mask.end(),
884 for (
unsigned int i = 1; i < n_lanes_filled; ++i)
892 typename Number::value_type,
893 Number>::apply(shape_info.n_components,
894 shape_info.data.front().fe_degree,
898 evaluation_data_coarse.
begin());
904 template <
int dim,
typename Number,
typename IndexType>
905 inline const typename Number::value_type *
907 const unsigned int row)
const
917 template <
int dim,
typename Number,
typename IndexType>
918 inline const typename Number::value_type *
920 const unsigned int row)
const
930 template <
int dim,
typename Number,
typename IndexType>
934 std::size_t size = 0;
963 template <
typename Number>
967 std::size_t size = 0;
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
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
ConstraintValues< double > constraint_values
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
IndexSet locally_owned_indices
std::vector< unsigned int > plain_dof_indices
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::size_t memory_consumption() const
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::vector< unsigned int > active_fe_indices
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
std::vector< unsigned int > dof_indices
void reinit(const unsigned int n_cells)
std::vector< ConstraintKinds > mask
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
#define DEAL_II_NAMESPACE_CLOSE
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)
std::uint8_t compressed_constraint_kind
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
static const unsigned int invalid_unsigned_int
const types::global_dof_index invalid_dof_index
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
std::size_t memory_consumption() const
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