Reference documentation for deal.II version Git e7bb9ce7b3 2020-09-18 12:07:32 -0400
\(\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 
19 
22 
23 #include <deal.II/fe/fe.h>
24 
25 #include <deal.II/grid/tria.h>
28 
32 #include <deal.II/lac/la_vector.h>
37 #include <deal.II/lac/vector.h>
39 
41 
43 
44 template <int dim, typename VectorType, typename DoFHandlerType>
46  const DoFHandlerType &dof)
47  : dof_handler(&dof, typeid(*this).name())
48  , n_dofs_old(0)
49  , prepared_for(none)
50 {
51  Assert((dynamic_cast<const parallel::distributed::Triangulation<
52  DoFHandlerType::dimension,
53  DoFHandlerType::space_dimension> *>(
54  &dof_handler->get_triangulation()) == nullptr),
55  ExcMessage("You are calling the ::SolutionTransfer class "
56  "with a DoF handler 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, typename DoFHandlerType>
66 {
67  clear();
68 }
69 
70 
71 
72 template <int dim, typename VectorType, typename DoFHandlerType>
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, typename DoFHandlerType>
86 void
88 {
92 
93  clear();
94 
95  const unsigned int n_active_cells =
96  dof_handler->get_triangulation().n_active_cells();
97  n_dofs_old = dof_handler->n_dofs();
98 
99  // efficient reallocation of indices_on_cell
100  std::vector<std::vector<types::global_dof_index>>(n_active_cells)
102 
103  typename DoFHandlerType::active_cell_iterator cell =
104  dof_handler->begin_active(),
105  endc = dof_handler->end();
106 
107  for (unsigned int i = 0; cell != endc; ++cell, ++i)
108  {
109  indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
110  // on each cell store the indices of the
111  // dofs. after refining we get the values
112  // on the children by taking these
113  // indices, getting the respective values
114  // out of the data vectors and prolonging
115  // them to the children
116  cell->get_dof_indices(indices_on_cell[i]);
117  cell_map[std::make_pair(cell->level(), cell->index())] =
118  Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
119  }
121 }
122 
123 
124 
125 template <int dim, typename VectorType, typename DoFHandlerType>
126 void
128  const VectorType &in,
129  VectorType & out) const
130 {
132  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
133  Assert(out.size() == dof_handler->n_dofs(),
134  ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
135  Assert(&in != &out,
136  ExcMessage("Vectors cannot be used as input and output"
137  " at the same time!"));
138 
140 
141  typename DoFHandlerType::cell_iterator cell = dof_handler->begin(),
142  endc = dof_handler->end();
143 
144  typename std::map<std::pair<unsigned int, unsigned int>,
145  Pointerstruct>::const_iterator pointerstruct,
146  cell_map_end = cell_map.end();
147 
148  for (; cell != endc; ++cell)
149  {
150  pointerstruct =
151  cell_map.find(std::make_pair(cell->level(), cell->index()));
152 
153  if (pointerstruct != cell_map_end)
154  // this cell was refined or not
155  // touched at all, so we can get
156  // the new values by just setting
157  // or interpolating to the children,
158  // which is both done by one
159  // function
160  {
161  const unsigned int this_fe_index =
162  pointerstruct->second.active_fe_index;
163  const unsigned int dofs_per_cell =
164  cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
165  local_values.reinit(dofs_per_cell, true);
166 
167  // make sure that the size of the stored indices is the same as
168  // dofs_per_cell. since we store the desired fe_index, we know
169  // what this size should be
170  Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
171  ExcInternalError());
172  for (unsigned int i = 0; i < dofs_per_cell; ++i)
174  in, (*pointerstruct->second.indices_ptr)[i]);
175  cell->set_dof_values_by_interpolation(local_values,
176  out,
177  this_fe_index);
178  }
179  }
180 }
181 
182 
183 
184 namespace internal
185 {
199  template <int dim, int spacedim>
200  void
201  extract_interpolation_matrices(const ::DoFHandler<dim, spacedim> &dof,
202  ::Table<2, FullMatrix<double>> &matrices)
203  {
204  if (dof.hp_capability_enabled == false)
205  return;
206 
207  const ::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
208  matrices.reinit(fe.size(), fe.size());
209  for (unsigned int i = 0; i < fe.size(); ++i)
210  for (unsigned int j = 0; j < fe.size(); ++j)
211  if (i != j)
212  {
213  matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
214  fe[j].n_dofs_per_cell());
215 
216  // see if we can get the interpolation matrices for this
217  // combination of elements. if not, reset the matrix sizes to zero
218  // to indicate that this particular combination isn't
219  // supported. this isn't an outright error right away since we may
220  // never need to actually interpolate between these two elements
221  // on actual cells; we simply have to trigger an error if someone
222  // actually tries
223  try
224  {
225  fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
226  }
227  catch (const typename FiniteElement<dim, spacedim>::
228  ExcInterpolationNotImplemented &)
229  {
230  matrices(i, j).reinit(0, 0);
231  }
232  }
233  }
234 
235 
236  template <int dim, int spacedim>
237  void
239  std::vector<std::vector<bool>> &)
240  {}
241 
242  template <int dim, int spacedim>
243  void
244  restriction_additive(const ::hp::FECollection<dim, spacedim> &fe,
245  std::vector<std::vector<bool>> &restriction_is_additive)
246  {
247  restriction_is_additive.resize(fe.size());
248  for (unsigned int f = 0; f < fe.size(); ++f)
249  {
250  restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
251  for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
252  restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
253  }
254  }
255 } // namespace internal
256 
257 
258 
259 template <int dim, typename VectorType, typename DoFHandlerType>
260 void
262  prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
263 {
267 
268  const unsigned int in_size = all_in.size();
269  Assert(in_size != 0,
270  ExcMessage("The array of input vectors you pass to this "
271  "function has no elements. This is not useful."));
272 
273  clear();
274 
275  const unsigned int n_active_cells =
276  dof_handler->get_triangulation().n_active_cells();
277  (void)n_active_cells;
278  n_dofs_old = dof_handler->n_dofs();
279 
280  for (unsigned int i = 0; i < in_size; ++i)
281  {
282  Assert(all_in[i].size() == n_dofs_old,
283  ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
284  }
285 
286  // first count the number
287  // of cells that will be coarsened
288  // and that'll stay or be refined
289  unsigned int n_cells_to_coarsen = 0;
290  unsigned int n_cells_to_stay_or_refine = 0;
291  for (const auto &act_cell : dof_handler->active_cell_iterators())
292  {
293  if (act_cell->coarsen_flag_set())
294  ++n_cells_to_coarsen;
295  else
296  ++n_cells_to_stay_or_refine;
297  }
298  Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) == n_active_cells,
299  ExcInternalError());
300 
301  unsigned int n_coarsen_fathers = 0;
302  for (typename DoFHandlerType::cell_iterator cell = dof_handler->begin();
303  cell != dof_handler->end();
304  ++cell)
305  if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
306  ++n_coarsen_fathers;
307  Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
308 
309  // allocate the needed memory. initialize
310  // the following arrays in an efficient
311  // way, without copying much
312  std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
314 
315  std::vector<std::vector<Vector<typename VectorType::value_type>>>(
316  n_coarsen_fathers,
317  std::vector<Vector<typename VectorType::value_type>>(in_size))
318  .swap(dof_values_on_cell);
319 
320  Table<2, FullMatrix<double>> interpolation_hp;
321  std::vector<std::vector<bool>> restriction_is_additive;
322 
324  internal::restriction_additive(dof_handler->get_fe_collection(),
325  restriction_is_additive);
326 
327  // we need counters for
328  // the 'to_stay_or_refine' cells 'n_sr' and
329  // the 'coarsen_fathers' cells 'n_cf',
330  unsigned int n_sr = 0, n_cf = 0;
331  for (typename DoFHandlerType::cell_iterator cell = dof_handler->begin();
332  cell != dof_handler->end();
333  ++cell)
334  {
335  // CASE 1: active cell that remains as it is
336  if (cell->is_active() && !cell->coarsen_flag_set())
337  {
338  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
339  indices_on_cell[n_sr].resize(dofs_per_cell);
340  // cell will not be coarsened,
341  // so we get away by storing the
342  // dof indices and later
343  // interpolating to the children
344  cell->get_dof_indices(indices_on_cell[n_sr]);
345  cell_map[std::make_pair(cell->level(), cell->index())] =
346  Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
347  ++n_sr;
348  }
349 
350  // CASE 2: cell is inactive but will become active
351  else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
352  {
353  // we will need to interpolate from the children of this cell
354  // to the current one. in the hp context, this also means
355  // we need to figure out which finite element space to interpolate
356  // to since that is not implied by the global FE as in the non-hp
357  // case. we choose the 'least dominant fe' on all children from
358  // the associated FECollection.
359  std::set<unsigned int> fe_indices_children;
360  for (unsigned int child_index = 0; child_index < cell->n_children();
361  ++child_index)
362  {
363  const auto &child = cell->child(child_index);
364  Assert(child->is_active() && child->coarsen_flag_set(),
366  dim>::ExcInconsistentCoarseningFlags());
367 
368  fe_indices_children.insert(child->active_fe_index());
369  }
370  Assert(!fe_indices_children.empty(), ExcInternalError());
371 
372  const unsigned int target_fe_index =
373  dof_handler->get_fe_collection().find_dominated_fe_extended(
374  fe_indices_children, /*codim=*/0);
375 
376  Assert(target_fe_index != numbers::invalid_unsigned_int,
377  typename ::hp::FECollection<
378  dim>::ExcNoDominatedFiniteElementAmongstChildren());
379 
380  const unsigned int dofs_per_cell =
381  dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
382 
383  std::vector<Vector<typename VectorType::value_type>>(
384  in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
385  .swap(dof_values_on_cell[n_cf]);
386 
387 
388  // store the data of each of the input vectors. get this data
389  // as interpolated onto a finite element space that encompasses
390  // that of all the children. note that
391  // cell->get_interpolated_dof_values already does all of the
392  // interpolations between spaces
393  for (unsigned int j = 0; j < in_size; ++j)
394  cell->get_interpolated_dof_values(all_in[j],
395  dof_values_on_cell[n_cf][j],
396  target_fe_index);
397  cell_map[std::make_pair(cell->level(), cell->index())] =
398  Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
399  ++n_cf;
400  }
401  }
402  Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
403  Assert(n_cf == n_coarsen_fathers, ExcInternalError());
404 
406 }
407 
408 
409 
410 template <int dim, typename VectorType, typename DoFHandlerType>
411 void
414 {
415  std::vector<VectorType> all_in(1, in);
417 }
418 
419 
420 
421 template <int dim, typename VectorType, typename DoFHandlerType>
422 void
424  const std::vector<VectorType> &all_in,
425  std::vector<VectorType> & all_out) const
426 {
428  const unsigned int size = all_in.size();
429  Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
430  for (unsigned int i = 0; i < size; ++i)
431  Assert(all_in[i].size() == n_dofs_old,
432  ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
433  for (unsigned int i = 0; i < all_out.size(); ++i)
434  Assert(all_out[i].size() == dof_handler->n_dofs(),
435  ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
436  for (unsigned int i = 0; i < size; ++i)
437  for (unsigned int j = 0; j < size; ++j)
438  Assert(&all_in[i] != &all_out[j],
439  ExcMessage("Vectors cannot be used as input and output"
440  " at the same time!"));
441 
443  std::vector<types::global_dof_index> dofs;
444 
445  typename std::map<std::pair<unsigned int, unsigned int>,
446  Pointerstruct>::const_iterator pointerstruct,
447  cell_map_end = cell_map.end();
448 
449  Table<2, FullMatrix<double>> interpolation_hp;
452 
453  for (const auto &cell : dof_handler->cell_iterators())
454  {
455  pointerstruct =
456  cell_map.find(std::make_pair(cell->level(), cell->index()));
457 
458  if (pointerstruct != cell_map_end)
459  {
460  const std::vector<types::global_dof_index> *const indexptr =
461  pointerstruct->second.indices_ptr;
462 
463  const std::vector<Vector<typename VectorType::value_type>>
464  *const valuesptr = pointerstruct->second.dof_values_ptr;
465 
466  // cell stayed as it was or was refined
467  if (indexptr)
468  {
469  Assert(valuesptr == nullptr, ExcInternalError());
470 
471  const unsigned int old_fe_index =
472  pointerstruct->second.active_fe_index;
473 
474  // get the values of each of the input data vectors on this cell
475  // and prolong it to its children
476  unsigned int in_size = indexptr->size();
477  for (unsigned int j = 0; j < size; ++j)
478  {
479  tmp.reinit(in_size, true);
480  for (unsigned int i = 0; i < in_size; ++i)
481  tmp(i) =
483  (*indexptr)[i]);
484 
485  cell->set_dof_values_by_interpolation(tmp,
486  all_out[j],
487  old_fe_index);
488  }
489  }
490  else if (valuesptr)
491  // the children of this cell were deleted
492  {
493  Assert(!cell->has_children(), ExcInternalError());
494  Assert(indexptr == nullptr, ExcInternalError());
495 
496  const unsigned int dofs_per_cell =
497  cell->get_fe().n_dofs_per_cell();
498  dofs.resize(dofs_per_cell);
499  // get the local
500  // indices
501  cell->get_dof_indices(dofs);
502 
503  // distribute the stored data to the new vectors
504  for (unsigned int j = 0; j < size; ++j)
505  {
506  // make sure that the size of the stored indices is the same
507  // as dofs_per_cell. this is kind of a test if we use the same
508  // fe in the hp case. to really do that test we would have to
509  // store the fe_index of all cells
510  const Vector<typename VectorType::value_type> *data = nullptr;
511  const unsigned int active_fe_index = cell->active_fe_index();
512  if (active_fe_index != pointerstruct->second.active_fe_index)
513  {
514  const unsigned int old_index =
515  pointerstruct->second.active_fe_index;
516  const FullMatrix<double> &interpolation_matrix =
517  interpolation_hp(active_fe_index, old_index);
518  // The interpolation matrix might be empty when using
519  // FE_Nothing.
520  if (interpolation_matrix.empty())
521  tmp.reinit(dofs_per_cell, false);
522  else
523  {
524  tmp.reinit(dofs_per_cell, true);
525  AssertDimension((*valuesptr)[j].size(),
526  interpolation_matrix.n());
527  AssertDimension(tmp.size(), interpolation_matrix.m());
528  interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
529  }
530  data = &tmp;
531  }
532  else
533  data = &(*valuesptr)[j];
534 
535 
536  for (unsigned int i = 0; i < dofs_per_cell; ++i)
538  dofs[i],
539  all_out[j]);
540  }
541  }
542  // undefined status
543  else
544  Assert(false, ExcInternalError());
545  }
546  }
547 }
548 
549 
550 
551 template <int dim, typename VectorType, typename DoFHandlerType>
552 void
554  const VectorType &in,
555  VectorType & out) const
556 {
557  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
558  Assert(out.size() == dof_handler->n_dofs(),
559  ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
560 
561  std::vector<VectorType> all_in(1);
562  all_in[0] = in;
563  std::vector<VectorType> all_out(1);
564  all_out[0] = out;
565  interpolate(all_in, all_out);
566  out = all_out[0];
567 }
568 
569 
570 
571 template <int dim, typename VectorType, typename DoFHandlerType>
572 std::size_t
574 {
575  // at the moment we do not include the memory
576  // consumption of the cell_map as we have no
577  // real idea about memory consumption of a
578  // std::map
581  sizeof(prepared_for) +
584 }
585 
586 
587 
588 template <int dim, typename VectorType, typename DoFHandlerType>
589 std::size_t
592 {
593  return sizeof(*this);
594 }
595 
596 
597 /*-------------- Explicit Instantiations -------------------------------*/
598 #define SPLIT_INSTANTIATIONS_COUNT 4
599 #ifndef SPLIT_INSTANTIATIONS_INDEX
600 # define SPLIT_INSTANTIATIONS_INDEX 0
601 #endif
602 #include "solution_transfer.inst"
603 
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
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:1411
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:363
void swap(SmartPointer< T, Q > &tt)
Definition: smartpointer.h:390
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:11751
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:302
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 > &)
Definition: memory_space.h:103
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:362
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
Definition: table.h:687
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()