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