Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.4.1
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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
31namespace 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.) <
253 100 * std::numeric_limits<double>::epsilon())
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,
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
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 hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
Definition: point.h:111
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
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
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:4604
unsigned int cell_index
Definition: grid_tools.cc:1129
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
compressed_constraint_kind compress(const ConstraintKinds kind_in, const unsigned int dim)
std::uint8_t compressed_constraint_kind
Definition: dof_info.h:75
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
static const unsigned int invalid_unsigned_int
Definition: types.h:201
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)