Reference documentation for deal.II version GIT 472adc65c1 2022-06-27 15:55:01+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\}}\)
constraint_info.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2022 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 
17 #ifndef dealii_matrix_free_constraint_info_h
18 #define dealii_matrix_free_constraint_info_h
19 
20 
21 #include <deal.II/base/config.h>
22 
27 
28 
30 
31 namespace internal
32 {
33  namespace MatrixFreeFunctions
34  {
40  template <int dim, typename Number>
42  {
43  public:
44  void
45  reinit(const DoFHandler<dim> &dof_handler,
46  const unsigned int n_cells,
47  const bool use_fast_hanging_node_algorithm = true);
48 
49  void
51  const unsigned int cell_no,
52  const unsigned int mg_level,
54  const ::AffineConstraints<typename Number::value_type>
55  & constraints,
56  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
57 
58  void
60 
61  template <typename T, typename VectorType>
62  void
63  read_write_operation(const T & operation,
64  VectorType & global_vector,
65  AlignedVector<Number> &local_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;
70 
71  void
73  const unsigned int first_cell,
74  const unsigned int n_lanes_filled,
75  const bool transpose,
76  AlignedVector<Number> &evaluation_data_coarse) const;
77 
78  private:
79  // for setup
81  std::vector<std::vector<unsigned int>> dof_indices_per_cell;
82  std::vector<std::vector<unsigned int>> plain_dof_indices_per_cell;
83  std::vector<std::vector<std::pair<unsigned short, unsigned short>>>
85 
86  std::unique_ptr<HangingNodes<dim>> hanging_nodes;
87  std::vector<std::vector<unsigned int>> lexicographic_numbering;
88 
89  public:
90  // for read_write_operation()
91  std::vector<unsigned int> dof_indices;
92  std::vector<std::pair<unsigned short, unsigned short>>
94  std::vector<std::pair<unsigned int, unsigned int>> row_starts;
95 
96  std::vector<unsigned int> plain_dof_indices;
97  std::vector<unsigned int> row_starts_plain_indices;
98 
99  // for constraint_pool_begin/end()
100  std::vector<typename Number::value_type> constraint_pool_data;
101  std::vector<unsigned int> constraint_pool_row_index;
102 
103  std::vector<ShapeInfo<Number>> shape_infos;
104  std::vector<compressed_constraint_kind> hanging_node_constraint_masks;
105  std::vector<unsigned int> active_fe_indices;
106 
107  private:
108  inline const typename Number::value_type *
109  constraint_pool_begin(const unsigned int row) const;
110 
111  inline const typename Number::value_type *
112  constraint_pool_end(const unsigned int row) const;
113  };
114 
115 
116 
117  template <int dim, typename Number>
118  inline void
120  const DoFHandler<dim> &dof_handler,
121  const unsigned int n_cells,
122  const bool use_fast_hanging_node_algorithm)
123  {
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);
127 
128  // note: has_hanging_nodes() is a global operatrion
129  const bool has_hanging_nodes =
130  dof_handler.get_triangulation().has_hanging_nodes();
131 
132  if (use_fast_hanging_node_algorithm && has_hanging_nodes)
133  {
134  hanging_nodes = std::make_unique<HangingNodes<dim>>(
135  dof_handler.get_triangulation());
136 
137  hanging_node_constraint_masks.resize(n_cells);
138  }
139 
140  const auto &fes = dof_handler.get_fe_collection();
141  lexicographic_numbering.resize(fes.size());
142  shape_infos.resize(fes.size());
143 
144  for (unsigned int i = 0; i < fes.size(); ++i)
145  {
146  if (fes[i].reference_cell().is_hyper_cube())
147  {
148  const Quadrature<1> dummy_quadrature(
149  std::vector<Point<1>>(1, Point<1>()));
150  shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
151  }
152  else
153  {
154  const auto dummy_quadrature =
155  fes[i].reference_cell().template get_gauss_type_quadrature<dim>(
156  1);
157  shape_infos[i].reinit(dummy_quadrature, fes[i], 0);
158  }
159 
160  lexicographic_numbering[i] = shape_infos[i].lexicographic_numbering;
161  }
162  active_fe_indices.resize(n_cells);
163  }
164 
165 
166 
167  template <int dim, typename Number>
168  inline void
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)
175  {
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());
180 
181  if (mg_level == numbers::invalid_unsigned_int)
182  cell->get_dof_indices(local_dof_indices);
183  else
184  cell->get_mg_dof_indices(local_dof_indices);
185 
186  {
187  AssertIndexRange(cell->active_fe_index(), shape_infos.size());
188 
189  const auto &lexicographic_numbering =
190  shape_infos[cell->active_fe_index()].lexicographic_numbering;
191 
192  AssertDimension(lexicographic_numbering.size(),
193  local_dof_indices.size());
194 
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]];
198  }
199 
200  std::pair<unsigned short, unsigned short> constraint_iterator(0, 0);
201 
202  AssertIndexRange(cell_no, this->constraint_indicator_per_cell.size());
203  AssertIndexRange(cell_no, this->dof_indices_per_cell.size());
204  AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size());
205  AssertIndexRange(cell_no, this->plain_dof_indices_per_cell.size());
206 
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];
210 
211  AssertDimension(constraint_indicator_per_cell[cell_no].size(), 0);
212  AssertDimension(dof_indices_per_cell[cell_no].size(), 0);
213  AssertDimension(plain_dof_indices_per_cell[cell_no].size(), 0);
214 
215  const auto global_to_local =
216  [&](const types::global_dof_index global_index) -> unsigned int {
217  if (partitioner)
218  return partitioner->global_to_local(global_index);
219  else
220  return global_index;
221  };
222 
223  // plain indices
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]);
227 
228  if (hanging_nodes)
229  {
230  AssertIndexRange(cell_no, this->hanging_node_constraint_masks.size());
231  AssertIndexRange(cell_no, this->active_fe_indices.size());
232 
233  std::vector<ConstraintKinds> mask(cell->get_fe().n_components());
234  hanging_nodes->setup_constraints(
235  cell, {}, lexicographic_numbering, local_dof_indices_lex, mask);
236 
237  hanging_node_constraint_masks[cell_no] = compress(mask[0], dim);
238  active_fe_indices[cell_no] = cell->active_fe_index();
239  }
240 
241  for (auto current_dof : local_dof_indices_lex)
242  {
243  const auto *entries_ptr =
244  constraints.get_constraint_entries(current_dof);
245 
246  // dof is constrained
247  if (entries_ptr != nullptr)
248  {
249  const auto & entries = *entries_ptr;
250  const types::global_dof_index n_entries = entries.size();
251  if (n_entries == 1 &&
252  std::abs(entries[0].second - 1.) <
254  {
255  current_dof = entries[0].first;
256  goto no_constraint;
257  }
258 
259  constraint_indicator.push_back(constraint_iterator);
260  constraint_indicator.back().second =
261  constraint_values.insert_entries(entries);
262 
263  // reset constraint iterator for next round
264  constraint_iterator.first = 0;
265 
266  if (n_entries > 0)
267  {
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)
271  {
272  dof_indices.push_back(
273  global_to_local(constraint_indices[j]));
274  }
275  }
276  }
277  else
278  {
279  no_constraint:
280  dof_indices.push_back(global_to_local(current_dof));
281 
282  // make sure constraint_iterator.first is always within the
283  // bounds of unsigned short
284  Assert(constraint_iterator.first <
285  (1 << (8 * sizeof(unsigned short))) - 1,
286  ExcInternalError());
287  constraint_iterator.first++;
288  }
289  }
290  }
291 
292 
293 
294  template <int dim, typename Number>
295  inline void
297  {
298  this->dof_indices = {};
299  this->plain_dof_indices = {};
300  this->constraint_indicator = {};
301 
302  this->row_starts = {};
303  this->row_starts.emplace_back(0, 0);
304 
305  this->row_starts_plain_indices = {};
306  this->row_starts_plain_indices.emplace_back(0);
307 
308  for (unsigned int i = 0; i < dof_indices_per_cell.size(); ++i)
309  {
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());
317 
318  this->row_starts.emplace_back(this->dof_indices.size(),
319  this->constraint_indicator.size());
320 
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());
324 
325  this->row_starts_plain_indices.emplace_back(
326  this->plain_dof_indices.size());
327  }
328 
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)
333  {
334  AssertIndexRange(it.second, constraints.size());
335  constraints[it.second] = &it.first;
336  length += it.first.size();
337  }
338 
339  constraint_pool_data.clear();
340  constraint_pool_data.reserve(length);
341  constraint_pool_row_index.reserve(constraint_values.constraints.size() +
342  1);
343  constraint_pool_row_index.resize(1, 0);
344 
345  for (const auto &constraint : constraints)
346  {
347  Assert(constraint != nullptr, ExcInternalError());
348  constraint_pool_data.insert(constraint_pool_data.end(),
349  constraint->begin(),
350  constraint->end());
351  constraint_pool_row_index.push_back(constraint_pool_data.size());
352  }
353 
354  AssertDimension(constraint_pool_data.size(), length);
355 
356  dof_indices_per_cell.clear();
357  constraint_indicator_per_cell.clear();
358 
359  if (hanging_nodes &&
360  std::all_of(hanging_node_constraint_masks.begin(),
361  hanging_node_constraint_masks.end(),
362  [](const auto i) {
363  return i == unconstrained_compressed_constraint_kind;
364  }))
365  hanging_node_constraint_masks.clear();
366  }
367 
368 
369 
370  template <int dim, typename Number>
371  template <typename T, typename VectorType>
372  inline void
374  const T & operation,
375  VectorType & global_vector,
376  AlignedVector<Number> &local_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
381  {
382  if (apply_constraints == false)
383  {
384  for (unsigned int v = 0; v < n_cells; ++v)
385  {
386  const unsigned int cell_index = first_cell + v;
387  const unsigned int *dof_indices =
388  this->plain_dof_indices.data() +
389  this->row_starts_plain_indices[cell_index];
390 
391  for (unsigned int i = 0; i < n_dofs_per_cell; ++dof_indices, ++i)
392  operation.process_dof(*dof_indices,
393  global_vector,
394  local_vector[i][v]);
395  }
396 
397  return;
398  }
399 
400  for (unsigned int v = 0; v < n_cells; ++v)
401  {
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 =
407  this->row_starts[cell_index + 1].second;
408 
409  unsigned int ind_local = 0;
410  for (; index_indicators != next_index_indicators; ++index_indicators)
411  {
412  const std::pair<unsigned short, unsigned short> indicator =
413  this->constraint_indicator[index_indicators];
414 
415  // run through values up to next constraint
416  for (unsigned int j = 0; j < indicator.first; ++j)
417  operation.process_dof(dof_indices[j],
418  global_vector,
419  local_vector[ind_local + j][v]);
420 
421  ind_local += indicator.first;
422  dof_indices += indicator.first;
423 
424  // constrained case: build the local value as a linear
425  // combination of the global value according to constraints
426  typename Number::value_type value;
427  operation.pre_constraints(local_vector[ind_local][v], value);
428 
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,
435  *data_val,
436  global_vector,
437  value);
438 
439  operation.post_constraints(value, local_vector[ind_local][v]);
440  ind_local++;
441  }
442 
443  AssertIndexRange(ind_local, n_dofs_per_cell + 1);
444 
445  for (; ind_local < n_dofs_per_cell; ++dof_indices, ++ind_local)
446  operation.process_dof(*dof_indices,
447  global_vector,
448  local_vector[ind_local][v]);
449  }
450  }
451 
452 
453 
454  template <int dim, typename Number>
455  inline void
457  const unsigned int first_cell,
458  const unsigned int n_lanes_filled,
459  const bool transpose,
460  AlignedVector<Number> &evaluation_data_coarse) const
461  {
462  if (hanging_node_constraint_masks.size() == 0)
463  return;
464 
466  Number::size()>
467  constraint_mask;
468 
469  bool hn_available = false;
470 
471  for (unsigned int v = 0; v < n_lanes_filled; ++v)
472  {
473  const auto mask = hanging_node_constraint_masks[first_cell + v];
474 
475  constraint_mask[v] = mask;
476 
477  hn_available |= (mask != unconstrained_compressed_constraint_kind);
478  }
479 
480  if (hn_available == true)
481  {
482  for (unsigned int v = n_lanes_filled; v < Number::size(); ++v)
483  constraint_mask[v] = unconstrained_compressed_constraint_kind;
484 
485  for (unsigned int i = 1; i < n_lanes_filled; ++i)
486  AssertDimension(active_fe_indices[first_cell],
487  active_fe_indices[first_cell + i]);
488 
489  const auto &shape_info = shape_infos[active_fe_indices[first_cell]];
490 
492  dim,
493  typename Number::value_type,
494  Number>::apply(shape_info.n_components,
495  shape_info.data.front().fe_degree,
496  shape_info,
497  transpose,
498  constraint_mask,
499  evaluation_data_coarse.begin());
500  }
501  }
502 
503 
504 
505  template <int dim, typename Number>
506  inline const typename Number::value_type *
508  const unsigned int row) const
509  {
510  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
511  return constraint_pool_data.empty() ?
512  nullptr :
513  constraint_pool_data.data() + constraint_pool_row_index[row];
514  }
515 
516 
517 
518  template <int dim, typename Number>
519  inline const typename Number::value_type *
521  const unsigned int row) const
522  {
523  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
524  return constraint_pool_data.empty() ?
525  nullptr :
526  constraint_pool_data.data() + constraint_pool_row_index[row + 1];
527  }
528 
529 
530  } // namespace MatrixFreeFunctions
531 } // namespace internal
532 
534 
535 #endif
iterator begin()
const Triangulation< dim, spacedim > & get_triangulation() const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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
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
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< 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
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
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition: grid_out.cc:4605
unsigned int cell_index
Definition: grid_tools.cc:1129
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char T
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
std::uint8_t compressed_constraint_kind
Definition: dof_info.h:83
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13734
static const unsigned int invalid_unsigned_int
Definition: types.h:201
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
unsigned int global_dof_index
Definition: types.h:76