Reference documentation for deal.II version Git 24fafe1087 2021-11-29 17:11:51 +0100
\(\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\}}\)
solution_transfer.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2021 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 
17 
20 
23 
24 #include <deal.II/fe/fe.h>
25 
26 #include <deal.II/grid/tria.h>
29 
33 #include <deal.II/lac/la_vector.h>
38 #include <deal.II/lac/vector.h>
40 
42 
44 
45 template <int dim, typename VectorType, int spacedim>
47  const DoFHandler<dim, spacedim> &dof)
48  : dof_handler(&dof, typeid(*this).name())
49  , n_dofs_old(0)
50  , prepared_for(none)
51 {
52  Assert(
54  &dof_handler->get_triangulation()) == nullptr),
55  ExcMessage("You are calling the ::SolutionTransfer class "
56  "with a DoFHandler that is built on a "
57  "parallel::distributed::Triangulation. This will not "
58  "work for parallel computations. You probably want to "
59  "use the parallel::distributed::SolutionTransfer class."));
60 }
61 
62 
63 
64 template <int dim, typename VectorType, int spacedim>
66 {
67  clear();
68 }
69 
70 
71 
72 template <int dim, typename VectorType, int spacedim>
73 void
75 {
76  indices_on_cell.clear();
77  dof_values_on_cell.clear();
78  cell_map.clear();
79 
81 }
82 
83 
84 
85 template <int dim, typename VectorType, int spacedim>
86 void
88 {
92 
93  clear();
94 
95  // We need to access dof indices on the entire domain. For
96  // parallel::shared::Triangulations, ownership of cells might change. If they
97  // allow artificial cells, we need to restore the "true" cell owners
98  // temporarily.
99  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
100  // the current set of subdomain ids, set subdomain ids to the "true" owner of
101  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
102  // and later restore these flags when it is destroyed.
104  spacedim>
105  subdomain_modifier(dof_handler->get_triangulation());
106 
107  const unsigned int n_active_cells =
108  dof_handler->get_triangulation().n_active_cells();
109  n_dofs_old = dof_handler->n_dofs();
110 
111  // efficient reallocation of indices_on_cell
112  std::vector<std::vector<types::global_dof_index>>(n_active_cells)
114 
115  for (const auto &cell : dof_handler->active_cell_iterators())
116  {
117  const unsigned int i = cell->active_cell_index();
118  indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
119  // on each cell store the indices of the
120  // dofs. after refining we get the values
121  // on the children by taking these
122  // indices, getting the respective values
123  // out of the data vectors and prolonging
124  // them to the children
125  cell->get_dof_indices(indices_on_cell[i]);
126  cell_map[std::make_pair(cell->level(), cell->index())] =
127  Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
128  }
130 }
131 
132 
133 
134 template <int dim, typename VectorType, int spacedim>
135 void
137  const VectorType &in,
138  VectorType & out) const
139 {
141  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
142  Assert(out.size() == dof_handler->n_dofs(),
143  ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
144  Assert(&in != &out,
145  ExcMessage("Vectors cannot be used as input and output"
146  " at the same time!"));
147 
148  // We need to access dof indices on the entire domain. For
149  // parallel::shared::Triangulations, ownership of cells might change. If they
150  // allow artificial cells, we need to restore the "true" cell owners
151  // temporarily.
152  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
153  // the current set of subdomain ids, set subdomain ids to the "true" owner of
154  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
155  // and later restore these flags when it is destroyed.
157  spacedim>
158  subdomain_modifier(dof_handler->get_triangulation());
159 
161 
162  typename std::map<std::pair<unsigned int, unsigned int>,
163  Pointerstruct>::const_iterator pointerstruct,
164  cell_map_end = cell_map.end();
165 
166  for (const auto &cell : dof_handler->cell_iterators())
167  {
168  pointerstruct =
169  cell_map.find(std::make_pair(cell->level(), cell->index()));
170 
171  if (pointerstruct != cell_map_end)
172  // this cell was refined or not
173  // touched at all, so we can get
174  // the new values by just setting
175  // or interpolating to the children,
176  // which is both done by one
177  // function
178  {
179  const unsigned int this_fe_index =
180  pointerstruct->second.active_fe_index;
181  const unsigned int dofs_per_cell =
182  cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
183  local_values.reinit(dofs_per_cell, true);
184 
185  // make sure that the size of the stored indices is the same as
186  // dofs_per_cell. since we store the desired fe_index, we know
187  // what this size should be
188  Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
189  ExcInternalError());
190  for (unsigned int i = 0; i < dofs_per_cell; ++i)
192  in, (*pointerstruct->second.indices_ptr)[i]);
193  cell->set_dof_values_by_interpolation(local_values,
194  out,
195  this_fe_index);
196  }
197  }
198 }
199 
200 
201 
202 namespace internal
203 {
217  template <int dim, int spacedim>
218  void
219  extract_interpolation_matrices(const ::DoFHandler<dim, spacedim> &dof,
220  ::Table<2, FullMatrix<double>> &matrices)
221  {
222  if (dof.has_hp_capabilities() == false)
223  return;
224 
225  const ::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
226  matrices.reinit(fe.size(), fe.size());
227  for (unsigned int i = 0; i < fe.size(); ++i)
228  for (unsigned int j = 0; j < fe.size(); ++j)
229  if (i != j)
230  {
231  matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
232  fe[j].n_dofs_per_cell());
233 
234  // see if we can get the interpolation matrices for this
235  // combination of elements. if not, reset the matrix sizes to zero
236  // to indicate that this particular combination isn't
237  // supported. this isn't an outright error right away since we may
238  // never need to actually interpolate between these two elements
239  // on actual cells; we simply have to trigger an error if someone
240  // actually tries
241  try
242  {
243  fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
244  }
245  catch (const typename FiniteElement<dim, spacedim>::
246  ExcInterpolationNotImplemented &)
247  {
248  matrices(i, j).reinit(0, 0);
249  }
250  }
251  }
252 
253 
254  template <int dim, int spacedim>
255  void
257  std::vector<std::vector<bool>> &)
258  {}
259 
260  template <int dim, int spacedim>
261  void
262  restriction_additive(const ::hp::FECollection<dim, spacedim> &fe,
263  std::vector<std::vector<bool>> &restriction_is_additive)
264  {
265  restriction_is_additive.resize(fe.size());
266  for (unsigned int f = 0; f < fe.size(); ++f)
267  {
268  restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
269  for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
270  restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
271  }
272  }
273 } // namespace internal
274 
275 
276 
277 template <int dim, typename VectorType, int spacedim>
278 void
280  prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
281 {
285 
286  clear();
287  n_dofs_old = dof_handler->n_dofs();
288  const unsigned int in_size = all_in.size();
289 
290 #ifdef DEBUG
291  Assert(in_size != 0,
292  ExcMessage("The array of input vectors you pass to this "
293  "function has no elements. This is not useful."));
294  for (unsigned int i = 0; i < in_size; ++i)
295  {
296  Assert(all_in[i].size() == n_dofs_old,
297  ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
298  }
299 #endif
300 
301  // We need to access dof indices on the entire domain. For
302  // parallel::shared::Triangulations, ownership of cells might change. If they
303  // allow artificial cells, we need to restore the "true" cell owners
304  // temporarily.
305  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
306  // the current set of subdomain ids, set subdomain ids to the "true" owner of
307  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
308  // and later restore these flags when it is destroyed.
310  spacedim>
311  subdomain_modifier(dof_handler->get_triangulation());
312 
313  // first count the number
314  // of cells that will be coarsened
315  // and that'll stay or be refined
316  unsigned int n_cells_to_coarsen = 0;
317  unsigned int n_cells_to_stay_or_refine = 0;
318  for (const auto &act_cell : dof_handler->active_cell_iterators())
319  {
320  if (act_cell->coarsen_flag_set())
321  ++n_cells_to_coarsen;
322  else
323  ++n_cells_to_stay_or_refine;
324  }
325  Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
326  dof_handler->get_triangulation().n_active_cells(),
327  ExcInternalError());
328 
329  unsigned int n_coarsen_fathers = 0;
330  for (const auto &cell : dof_handler->cell_iterators())
331  if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
332  ++n_coarsen_fathers;
333  Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
334 
335  // allocate the needed memory. initialize
336  // the following arrays in an efficient
337  // way, without copying much
338  std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
340 
341  std::vector<std::vector<Vector<typename VectorType::value_type>>>(
342  n_coarsen_fathers,
343  std::vector<Vector<typename VectorType::value_type>>(in_size))
344  .swap(dof_values_on_cell);
345 
346  Table<2, FullMatrix<double>> interpolation_hp;
347  std::vector<std::vector<bool>> restriction_is_additive;
348 
350  internal::restriction_additive(dof_handler->get_fe_collection(),
351  restriction_is_additive);
352 
353  // we need counters for
354  // the 'to_stay_or_refine' cells 'n_sr' and
355  // the 'coarsen_fathers' cells 'n_cf',
356  unsigned int n_sr = 0, n_cf = 0;
357  for (const auto &cell : dof_handler->cell_iterators())
358  {
359  // CASE 1: active cell that remains as it is
360  if (cell->is_active() && !cell->coarsen_flag_set())
361  {
362  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
363  indices_on_cell[n_sr].resize(dofs_per_cell);
364  // cell will not be coarsened,
365  // so we get away by storing the
366  // dof indices and later
367  // interpolating to the children
368  cell->get_dof_indices(indices_on_cell[n_sr]);
369  cell_map[std::make_pair(cell->level(), cell->index())] =
370  Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
371  ++n_sr;
372  }
373 
374  // CASE 2: cell is inactive but will become active
375  else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
376  {
377  // we will need to interpolate from the children of this cell
378  // to the current one. in the hp-context, this also means
379  // we need to figure out which finite element space to interpolate
380  // to since that is not implied by the global FE as in the non-hp-
381  // case. we choose the 'least dominant fe' on all children from
382  // the associated FECollection.
383  std::set<unsigned int> fe_indices_children;
384  for (const auto &child : cell->child_iterators())
385  {
386  Assert(child->is_active() && child->coarsen_flag_set(),
387  typename ::Triangulation<
388  dim>::ExcInconsistentCoarseningFlags());
389 
390  fe_indices_children.insert(child->active_fe_index());
391  }
392  Assert(!fe_indices_children.empty(), ExcInternalError());
393 
394  const unsigned int target_fe_index =
395  dof_handler->get_fe_collection().find_dominated_fe_extended(
396  fe_indices_children, /*codim=*/0);
397 
398  Assert(target_fe_index != numbers::invalid_unsigned_int,
401 
402  const unsigned int dofs_per_cell =
403  dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
404 
405  std::vector<Vector<typename VectorType::value_type>>(
406  in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
407  .swap(dof_values_on_cell[n_cf]);
408 
409 
410  // store the data of each of the input vectors. get this data
411  // as interpolated onto a finite element space that encompasses
412  // that of all the children. note that
413  // cell->get_interpolated_dof_values already does all of the
414  // interpolations between spaces
415  for (unsigned int j = 0; j < in_size; ++j)
416  cell->get_interpolated_dof_values(all_in[j],
417  dof_values_on_cell[n_cf][j],
418  target_fe_index);
419  cell_map[std::make_pair(cell->level(), cell->index())] =
420  Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
421  ++n_cf;
422  }
423  }
424  Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
425  Assert(n_cf == n_coarsen_fathers, ExcInternalError());
426 
428 }
429 
430 
431 
432 template <int dim, typename VectorType, int spacedim>
433 void
436 {
437  std::vector<VectorType> all_in(1, in);
439 }
440 
441 
442 
443 template <int dim, typename VectorType, int spacedim>
444 void
446  const std::vector<VectorType> &all_in,
447  std::vector<VectorType> & all_out) const
448 {
449  const unsigned int size = all_in.size();
450 #ifdef DEBUG
452  Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
453  for (unsigned int i = 0; i < size; ++i)
454  Assert(all_in[i].size() == n_dofs_old,
455  ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
456  for (unsigned int i = 0; i < all_out.size(); ++i)
457  Assert(all_out[i].size() == dof_handler->n_dofs(),
458  ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
459  for (unsigned int i = 0; i < size; ++i)
460  for (unsigned int j = 0; j < size; ++j)
461  Assert(&all_in[i] != &all_out[j],
462  ExcMessage("Vectors cannot be used as input and output"
463  " at the same time!"));
464 #endif
465 
466  // We need to access dof indices on the entire domain. For
467  // parallel::shared::Triangulations, ownership of cells might change. If they
468  // allow artificial cells, we need to restore the "true" cell owners
469  // temporarily.
470  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
471  // the current set of subdomain ids, set subdomain ids to the "true" owner of
472  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
473  // and later restore these flags when it is destroyed.
475  spacedim>
476  subdomain_modifier(dof_handler->get_triangulation());
477 
479  std::vector<types::global_dof_index> dofs;
480 
481  typename std::map<std::pair<unsigned int, unsigned int>,
482  Pointerstruct>::const_iterator pointerstruct,
483  cell_map_end = cell_map.end();
484 
485  Table<2, FullMatrix<double>> interpolation_hp;
488 
489  for (const auto &cell : dof_handler->cell_iterators())
490  {
491  pointerstruct =
492  cell_map.find(std::make_pair(cell->level(), cell->index()));
493 
494  if (pointerstruct != cell_map_end)
495  {
496  const std::vector<types::global_dof_index> *const indexptr =
497  pointerstruct->second.indices_ptr;
498 
499  const std::vector<Vector<typename VectorType::value_type>>
500  *const valuesptr = pointerstruct->second.dof_values_ptr;
501 
502  // cell stayed as it was or was refined
503  if (indexptr)
504  {
505  Assert(valuesptr == nullptr, ExcInternalError());
506 
507  const unsigned int old_fe_index =
508  pointerstruct->second.active_fe_index;
509 
510  // get the values of each of the input data vectors on this cell
511  // and prolong it to its children
512  unsigned int in_size = indexptr->size();
513  for (unsigned int j = 0; j < size; ++j)
514  {
515  tmp.reinit(in_size, true);
516  for (unsigned int i = 0; i < in_size; ++i)
517  tmp(i) =
519  (*indexptr)[i]);
520 
521  cell->set_dof_values_by_interpolation(tmp,
522  all_out[j],
523  old_fe_index);
524  }
525  }
526  else if (valuesptr)
527  // the children of this cell were deleted
528  {
529  Assert(!cell->has_children(), ExcInternalError());
530  Assert(indexptr == nullptr, ExcInternalError());
531 
532  const unsigned int dofs_per_cell =
533  cell->get_fe().n_dofs_per_cell();
534  dofs.resize(dofs_per_cell);
535  // get the local
536  // indices
537  cell->get_dof_indices(dofs);
538 
539  // distribute the stored data to the new vectors
540  for (unsigned int j = 0; j < size; ++j)
541  {
542  // make sure that the size of the stored indices is the same
543  // as dofs_per_cell. this is kind of a test if we use the same
544  // FE in the hp-case. to really do that test we would have to
545  // store the fe_index of all cells
546  const Vector<typename VectorType::value_type> *data = nullptr;
547  const unsigned int active_fe_index = cell->active_fe_index();
548  if (active_fe_index != pointerstruct->second.active_fe_index)
549  {
550  const unsigned int old_index =
551  pointerstruct->second.active_fe_index;
552  const FullMatrix<double> &interpolation_matrix =
553  interpolation_hp(active_fe_index, old_index);
554  // The interpolation matrix might be empty when using
555  // FE_Nothing.
556  if (interpolation_matrix.empty())
557  tmp.reinit(dofs_per_cell, false);
558  else
559  {
560  tmp.reinit(dofs_per_cell, true);
561  AssertDimension((*valuesptr)[j].size(),
562  interpolation_matrix.n());
563  AssertDimension(tmp.size(), interpolation_matrix.m());
564  interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
565  }
566  data = &tmp;
567  }
568  else
569  data = &(*valuesptr)[j];
570 
571 
572  for (unsigned int i = 0; i < dofs_per_cell; ++i)
574  dofs[i],
575  all_out[j]);
576  }
577  }
578  // undefined status
579  else
580  Assert(false, ExcInternalError());
581  }
582  }
583 }
584 
585 
586 
587 template <int dim, typename VectorType, int spacedim>
588 void
590  VectorType &out) const
591 {
592  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
593  Assert(out.size() == dof_handler->n_dofs(),
594  ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
595 
596  std::vector<VectorType> all_in(1);
597  all_in[0] = in;
598  std::vector<VectorType> all_out(1);
599  all_out[0] = out;
600  interpolate(all_in, all_out);
601  out = all_out[0];
602 }
603 
604 
605 
606 template <int dim, typename VectorType, int spacedim>
607 std::size_t
609 {
610  // at the moment we do not include the memory
611  // consumption of the cell_map as we have no
612  // real idea about memory consumption of a
613  // std::map
616  sizeof(prepared_for) +
619 }
620 
621 
622 
623 template <int dim, typename VectorType, int spacedim>
624 std::size_t
626  const
627 {
628  return sizeof(*this);
629 }
630 
631 
632 /*-------------- Explicit Instantiations -------------------------------*/
633 #define SPLIT_INSTANTIATIONS_COUNT 4
634 #ifndef SPLIT_INSTANTIATIONS_INDEX
635 # define SPLIT_INSTANTIATIONS_INDEX 0
636 #endif
637 #include "solution_transfer.inst"
638 
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
static ::ExceptionBase & ExcNoDominatedFiniteElementOnChildren()
void prepare_for_pure_refinement()
PreparationState prepared_for
void interpolate(const std::vector< VectorType > &all_in, std::vector< VectorType > &all_out) const
static ::ExceptionBase & ExcNotPrepared()
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
std::map< std::pair< unsigned int, unsigned int >, Pointerstruct > cell_map
void extract_interpolation_matrices(const ::DoFHandler< dim, spacedim > &dof, ::Table< 2, FullMatrix< double >> &matrices)
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void restriction_additive(const FiniteElement< dim, spacedim > &, std::vector< std::vector< bool >> &)
SolutionTransfer(const DoFHandler< dim, spacedim > &dof)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
void swap(SmartPointer< T, Q > &tt)
Definition: smartpointer.h:394
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13728
std::vector< std::vector< types::global_dof_index > > indices_on_cell
std::size_t memory_consumption() const
void refine_interpolate(const VectorType &in, VectorType &out) const
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
SmartPointer< const DoFHandler< dim, spacedim >, SolutionTransfer< dim, VectorType, spacedim > > dof_handler
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
static ::ExceptionBase & ExcAlreadyPrepForRef()
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
std::size_t memory_consumption() const
std::vector< std::vector< Vector< typename VectorType::value_type > > > dof_values_on_cell
void prepare_for_coarsening_and_refinement(const std::vector< VectorType > &all_in)
types::global_dof_index n_dofs_old
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
void restriction_additive(const ::hp::FECollection< dim, spacedim > &fe, std::vector< std::vector< bool >> &restriction_is_additive)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAlreadyPrepForCoarseAndRef()