17#ifndef dealii_matrix_free_constraint_info_h
18#define dealii_matrix_free_constraint_info_h
33 namespace MatrixFreeFunctions
40 template <
int dim,
typename Number>
46 const unsigned int n_cells,
47 const bool use_fast_hanging_node_algorithm =
true);
51 const unsigned int cell_no,
52 const unsigned int mg_level,
54 const ::AffineConstraints<typename Number::value_type>
56 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
61 template <
typename T,
typename VectorType>
64 VectorType & global_vector,
66 const unsigned int first_cell,
67 const unsigned int n_cells,
68 const unsigned int n_dofs_per_cell,
69 const bool apply_constraints)
const;
73 const unsigned int first_cell,
74 const unsigned int n_lanes_filled,
83 std::vector<std::vector<std::pair<unsigned short, unsigned short>>>
92 std::vector<std::pair<unsigned short, unsigned short>>
94 std::vector<std::pair<unsigned int, unsigned int>>
row_starts;
108 inline const typename Number::value_type *
111 inline const typename Number::value_type *
117 template <
int dim,
typename Number>
121 const unsigned int n_cells,
122 const bool use_fast_hanging_node_algorithm)
124 this->dof_indices_per_cell.resize(n_cells);
125 this->plain_dof_indices_per_cell.resize(n_cells);
126 this->constraint_indicator_per_cell.resize(n_cells);
129 const bool has_hanging_nodes =
132 if (use_fast_hanging_node_algorithm && has_hanging_nodes)
134 hanging_nodes = std::make_unique<HangingNodes<dim>>(
137 hanging_node_constraint_masks.resize(n_cells);
141 lexicographic_numbering.resize(fes.size());
142 shape_infos.resize(fes.size());
144 for (
unsigned int i = 0; i < fes.size(); ++i)
146 if (fes[i].reference_cell().is_hyper_cube())
150 shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
154 const auto dummy_quadrature =
155 fes[i].reference_cell().template get_gauss_type_quadrature<dim>(
157 shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
160 lexicographic_numbering[i] = shape_infos[i].lexicographic_numbering;
162 active_fe_indices.resize(n_cells);
167 template <
int dim,
typename Number>
170 const unsigned int cell_no,
171 const unsigned int mg_level,
173 const ::AffineConstraints<typename Number::value_type> &constraints,
174 const std::shared_ptr<const Utilities::MPI::Partitioner> & partitioner)
176 std::vector<types::global_dof_index> local_dof_indices(
177 cell->get_fe().n_dofs_per_cell());
178 std::vector<types::global_dof_index> local_dof_indices_lex(
179 cell->get_fe().n_dofs_per_cell());
182 cell->get_dof_indices(local_dof_indices);
184 cell->get_mg_dof_indices(local_dof_indices);
189 const auto &lexicographic_numbering =
190 shape_infos[cell->active_fe_index()].lexicographic_numbering;
193 local_dof_indices.size());
195 for (
unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
196 local_dof_indices_lex[i] =
197 local_dof_indices[lexicographic_numbering[i]];
200 std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
207 auto &constraint_indicator = this->constraint_indicator_per_cell[cell_no];
208 auto &dof_indices = this->dof_indices_per_cell[cell_no];
209 auto &plain_dof_indices = this->plain_dof_indices_per_cell[cell_no];
215 const auto global_to_local =
218 return partitioner->global_to_local(global_index);
224 plain_dof_indices.resize(local_dof_indices_lex.size());
225 for (
unsigned int i = 0; i < local_dof_indices_lex.size(); ++i)
226 plain_dof_indices[i] = global_to_local(local_dof_indices_lex[i]);
233 std::vector<ConstraintKinds>
mask(cell->get_fe().n_components());
234 hanging_nodes->setup_constraints(
235 cell, {}, lexicographic_numbering, local_dof_indices_lex,
mask);
237 hanging_node_constraint_masks[cell_no] =
compress(
mask[0], dim);
238 active_fe_indices[cell_no] = cell->active_fe_index();
241 for (
auto current_dof : local_dof_indices_lex)
243 const auto *entries_ptr =
244 constraints.get_constraint_entries(current_dof);
247 if (entries_ptr !=
nullptr)
249 const auto & entries = *entries_ptr;
251 if (n_entries == 1 &&
253 100 * std::numeric_limits<double>::epsilon())
255 current_dof = entries[0].first;
259 constraint_indicator.push_back(constraint_iterator);
260 constraint_indicator.back().second =
261 constraint_values.insert_entries(entries);
264 constraint_iterator.first = 0;
268 const std::vector<types::global_dof_index>
269 &constraint_indices = constraint_values.constraint_indices;
270 for (
unsigned int j = 0; j < n_entries; ++j)
272 dof_indices.push_back(
273 global_to_local(constraint_indices[j]));
280 dof_indices.push_back(global_to_local(current_dof));
284 Assert(constraint_iterator.first <
285 (1 << (8 *
sizeof(
unsigned short))) - 1,
287 constraint_iterator.first++;
294 template <
int dim,
typename Number>
298 this->dof_indices = {};
299 this->plain_dof_indices = {};
300 this->constraint_indicator = {};
302 this->row_starts = {};
303 this->row_starts.emplace_back(0, 0);
305 this->row_starts_plain_indices = {};
306 this->row_starts_plain_indices.emplace_back(0);
308 for (
unsigned int i = 0; i < dof_indices_per_cell.size(); ++i)
310 this->dof_indices.insert(this->dof_indices.end(),
311 dof_indices_per_cell[i].begin(),
312 dof_indices_per_cell[i].end());
313 this->constraint_indicator.insert(
314 this->constraint_indicator.end(),
315 constraint_indicator_per_cell[i].begin(),
316 constraint_indicator_per_cell[i].end());
318 this->row_starts.emplace_back(this->dof_indices.size(),
319 this->constraint_indicator.size());
321 this->plain_dof_indices.insert(this->plain_dof_indices.end(),
322 plain_dof_indices_per_cell[i].begin(),
323 plain_dof_indices_per_cell[i].end());
325 this->row_starts_plain_indices.emplace_back(
326 this->plain_dof_indices.size());
329 std::vector<const std::vector<double> *> constraints(
330 constraint_values.constraints.size());
331 unsigned int length = 0;
332 for (
const auto &it : constraint_values.constraints)
335 constraints[it.second] = &it.first;
336 length += it.first.size();
339 constraint_pool_data.clear();
340 constraint_pool_data.reserve(length);
341 constraint_pool_row_index.reserve(constraint_values.constraints.size() +
343 constraint_pool_row_index.resize(1, 0);
345 for (
const auto &constraint : constraints)
348 constraint_pool_data.insert(constraint_pool_data.end(),
351 constraint_pool_row_index.push_back(constraint_pool_data.size());
356 dof_indices_per_cell.clear();
357 constraint_indicator_per_cell.clear();
360 std::all_of(hanging_node_constraint_masks.begin(),
361 hanging_node_constraint_masks.end(),
363 return i == unconstrained_compressed_constraint_kind;
365 hanging_node_constraint_masks.clear();
370 template <
int dim,
typename Number>
371 template <
typename T,
typename VectorType>
375 VectorType & global_vector,
377 const unsigned int first_cell,
378 const unsigned int n_cells,
379 const unsigned int n_dofs_per_cell,
380 const bool apply_constraints)
const
382 if (apply_constraints ==
false)
384 for (
unsigned int v = 0; v < n_cells; ++v)
386 const unsigned int cell_index = first_cell + v;
387 const unsigned int *dof_indices =
388 this->plain_dof_indices.data() +
391 for (
unsigned int i = 0; i < n_dofs_per_cell; ++dof_indices, ++i)
392 operation.process_dof(*dof_indices,
400 for (
unsigned int v = 0; v < n_cells; ++v)
402 const unsigned int cell_index = first_cell + v;
403 const unsigned int *dof_indices =
404 this->dof_indices.data() + this->row_starts[
cell_index].first;
405 unsigned int index_indicators = this->row_starts[
cell_index].second;
406 unsigned int next_index_indicators =
409 unsigned int ind_local = 0;
410 for (; index_indicators != next_index_indicators; ++index_indicators)
412 const std::pair<unsigned short, unsigned short> indicator =
413 this->constraint_indicator[index_indicators];
416 for (
unsigned int j = 0; j < indicator.first; ++j)
417 operation.process_dof(dof_indices[j],
419 local_vector[ind_local + j][v]);
421 ind_local += indicator.first;
422 dof_indices += indicator.first;
426 typename Number::value_type
value;
427 operation.pre_constraints(local_vector[ind_local][v],
value);
429 const typename Number::value_type *data_val =
430 this->constraint_pool_begin(indicator.second);
431 const typename Number::value_type *end_pool =
432 this->constraint_pool_end(indicator.second);
433 for (; data_val != end_pool; ++data_val, ++dof_indices)
434 operation.process_constraint(*dof_indices,
439 operation.post_constraints(
value, local_vector[ind_local][v]);
445 for (; ind_local < n_dofs_per_cell; ++dof_indices, ++ind_local)
446 operation.process_dof(*dof_indices,
448 local_vector[ind_local][v]);
454 template <
int dim,
typename Number>
457 const unsigned int first_cell,
458 const unsigned int n_lanes_filled,
462 if (hanging_node_constraint_masks.size() == 0)
469 bool hn_available =
false;
471 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
473 const auto mask = hanging_node_constraint_masks[first_cell + v];
475 constraint_mask[v] =
mask;
480 if (hn_available ==
true)
482 for (
unsigned int v = n_lanes_filled; v < Number::size(); ++v)
485 for (
unsigned int i = 1; i < n_lanes_filled; ++i)
487 active_fe_indices[first_cell + i]);
489 const auto &shape_info = shape_infos[active_fe_indices[first_cell]];
493 typename Number::value_type,
494 Number>::apply(shape_info.n_components,
495 shape_info.data.front().fe_degree,
499 evaluation_data_coarse.
begin());
505 template <
int dim,
typename Number>
506 inline const typename Number::value_type *
508 const unsigned int row)
const
511 return constraint_pool_data.empty() ?
513 constraint_pool_data.data() + constraint_pool_row_index[row];
518 template <
int dim,
typename Number>
519 inline const typename Number::value_type *
521 const unsigned int row)
const
524 return constraint_pool_data.empty() ?
526 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
virtual bool has_hanging_nodes() const
std::vector< std::vector< unsigned int > > lexicographic_numbering
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
void reinit(const DoFHandler< dim > &dof_handler, const unsigned int n_cells, const bool use_fast_hanging_node_algorithm=true)
std::vector< unsigned int > plain_dof_indices
std::vector< std::vector< std::pair< unsigned short, unsigned short > > > constraint_indicator_per_cell
std::vector< unsigned int > constraint_pool_row_index
std::vector< unsigned int > dof_indices
const Number::value_type * constraint_pool_begin(const unsigned int row) const
std::vector< unsigned int > active_fe_indices
std::vector< std::vector< unsigned int > > plain_dof_indices_per_cell
std::vector< std::vector< unsigned int > > dof_indices_per_cell
std::vector< ShapeInfo< Number > > shape_infos
ConstraintValues< double > constraint_values
std::unique_ptr< HangingNodes< dim > > hanging_nodes
std::vector< unsigned int > row_starts_plain_indices
const Number::value_type * constraint_pool_end(const unsigned int row) const
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
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
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::pair< unsigned int, unsigned int > > row_starts
std::vector< typename Number::value_type > constraint_pool_data
void read_write_operation(const T &operation, VectorType &global_vector, AlignedVector< 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
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
std::uint8_t compressed_constraint_kind
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)