Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solution_transfer.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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 #include <deal.II/base/memory_consumption.h>
17 
18 #include <deal.II/distributed/tria.h>
19 
20 #include <deal.II/dofs/dof_accessor.h>
21 #include <deal.II/dofs/dof_handler.h>
22 
23 #include <deal.II/fe/fe.h>
24 
25 #include <deal.II/grid/tria.h>
26 #include <deal.II/grid/tria_accessor.h>
27 #include <deal.II/grid/tria_iterator.h>
28 
29 #include <deal.II/lac/block_vector.h>
30 #include <deal.II/lac/la_parallel_block_vector.h>
31 #include <deal.II/lac/la_parallel_vector.h>
32 #include <deal.II/lac/la_vector.h>
33 #include <deal.II/lac/petsc_block_vector.h>
34 #include <deal.II/lac/petsc_vector.h>
35 #include <deal.II/lac/trilinos_parallel_block_vector.h>
36 #include <deal.II/lac/trilinos_vector.h>
37 #include <deal.II/lac/vector.h>
38 #include <deal.II/lac/vector_element_access.h>
39 
40 #include <deal.II/numerics/solution_transfer.h>
41 
42 DEAL_II_NAMESPACE_OPEN
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().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).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)
173  local_values(i) = internal::ElementAccess<VectorType>::get(
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 <typename DoFHandlerType>
200  void
201  extract_interpolation_matrices(const DoFHandlerType &,
202  ::Table<2, FullMatrix<double>> &)
203  {}
204 
205  template <int dim, int spacedim>
206  void
208  const ::hp::DoFHandler<dim, spacedim> &dof,
209  ::Table<2, FullMatrix<double>> & matrices)
210  {
211  const ::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
212  matrices.reinit(fe.size(), fe.size());
213  for (unsigned int i = 0; i < fe.size(); ++i)
214  for (unsigned int j = 0; j < fe.size(); ++j)
215  if (i != j)
216  {
217  matrices(i, j).reinit(fe[i].dofs_per_cell, fe[j].dofs_per_cell);
218 
219  // see if we can get the interpolation matrices for this
220  // combination of elements. if not, reset the matrix sizes to zero
221  // to indicate that this particular combination isn't
222  // supported. this isn't an outright error right away since we may
223  // never need to actually interpolate between these two elements
224  // on actual cells; we simply have to trigger an error if someone
225  // actually tries
226  try
227  {
228  fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
229  }
230  catch (const typename FiniteElement<dim, spacedim>::
231  ExcInterpolationNotImplemented &)
232  {
233  matrices(i, j).reinit(0, 0);
234  }
235  }
236  }
237 
238 
239  template <int dim, int spacedim>
240  void
241  restriction_additive(const FiniteElement<dim, spacedim> &,
242  std::vector<std::vector<bool>> &)
243  {}
244 
245  template <int dim, int spacedim>
246  void
247  restriction_additive(const ::hp::FECollection<dim, spacedim> &fe,
248  std::vector<std::vector<bool>> &restriction_is_additive)
249  {
250  restriction_is_additive.resize(fe.size());
251  for (unsigned int f = 0; f < fe.size(); ++f)
252  {
253  restriction_is_additive[f].resize(fe[f].dofs_per_cell);
254  for (unsigned int i = 0; i < fe[f].dofs_per_cell; ++i)
255  restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
256  }
257  }
258 } // namespace internal
259 
260 
261 
262 template <int dim, typename VectorType, typename DoFHandlerType>
263 void
265  prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
266 {
270 
271  const unsigned int in_size = all_in.size();
272  Assert(in_size != 0,
273  ExcMessage("The array of input vectors you pass to this "
274  "function has no elements. This is not useful."));
275 
276  clear();
277 
278  const unsigned int n_active_cells =
279  dof_handler->get_triangulation().n_active_cells();
280  (void)n_active_cells;
281  n_dofs_old = dof_handler->n_dofs();
282 
283  for (unsigned int i = 0; i < in_size; ++i)
284  {
285  Assert(all_in[i].size() == n_dofs_old,
286  ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
287  }
288 
289  // first count the number
290  // of cells that will be coarsened
291  // and that'll stay or be refined
292  unsigned int n_cells_to_coarsen = 0;
293  unsigned int n_cells_to_stay_or_refine = 0;
294  for (typename DoFHandlerType::active_cell_iterator act_cell =
295  dof_handler->begin_active();
296  act_cell != dof_handler->end();
297  ++act_cell)
298  {
299  if (act_cell->coarsen_flag_set())
300  ++n_cells_to_coarsen;
301  else
302  ++n_cells_to_stay_or_refine;
303  }
304  Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) == n_active_cells,
305  ExcInternalError());
306 
307  unsigned int n_coarsen_fathers = 0;
308  for (typename DoFHandlerType::cell_iterator cell = dof_handler->begin();
309  cell != dof_handler->end();
310  ++cell)
311  if (!cell->active() && cell->child(0)->coarsen_flag_set())
312  ++n_coarsen_fathers;
313  Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
314 
315  // allocate the needed memory. initialize
316  // the following arrays in an efficient
317  // way, without copying much
318  std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
320 
321  std::vector<std::vector<Vector<typename VectorType::value_type>>>(
322  n_coarsen_fathers,
323  std::vector<Vector<typename VectorType::value_type>>(in_size))
324  .swap(dof_values_on_cell);
325 
326  Table<2, FullMatrix<double>> interpolation_hp;
327  std::vector<std::vector<bool>> restriction_is_additive;
328 
330  internal::restriction_additive(dof_handler->get_fe_collection(),
331  restriction_is_additive);
332 
333  // we need counters for
334  // the 'to_stay_or_refine' cells 'n_sr' and
335  // the 'coarsen_fathers' cells 'n_cf',
336  unsigned int n_sr = 0, n_cf = 0;
337  for (typename DoFHandlerType::cell_iterator cell = dof_handler->begin();
338  cell != dof_handler->end();
339  ++cell)
340  {
341  // CASE 1: active cell that remains as it is
342  if (cell->active() && !cell->coarsen_flag_set())
343  {
344  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
345  indices_on_cell[n_sr].resize(dofs_per_cell);
346  // cell will not be coarsened,
347  // so we get away by storing the
348  // dof indices and later
349  // interpolating to the children
350  cell->get_dof_indices(indices_on_cell[n_sr]);
351  cell_map[std::make_pair(cell->level(), cell->index())] =
352  Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
353  ++n_sr;
354  }
355 
356  // CASE 2: cell is inactive but will become active
357  else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
358  {
359  // note that if one child has the
360  // coarsen flag, then all should
361  // have if Tria::prepare_* has
362  // worked correctly
363  for (unsigned int i = 1; i < cell->n_children(); ++i)
364  Assert(cell->child(i)->coarsen_flag_set(),
365  ExcMessage(
366  "It looks like you didn't call "
367  "Triangulation::prepare_coarsening_and_refinement before "
368  "calling the current function. This can't work."));
369 
370  // we will need to interpolate from the children of this cell
371  // to the current one. in the hp context, this also means
372  // we need to figure out which finite element space to interpolate
373  // to since that is not implied by the global FE as in the non-hp
374  // case.
375  bool different_fe_on_children = false;
376  for (unsigned int child = 1; child < cell->n_children(); ++child)
377  if (cell->child(child)->active_fe_index() !=
378  cell->child(0)->active_fe_index())
379  {
380  different_fe_on_children = true;
381  break;
382  }
383 
384  // take FE index from the child with most
385  // degrees of freedom locally
386  unsigned int most_general_child = 0;
387  if (different_fe_on_children == true)
388  for (unsigned int child = 1; child < cell->n_children(); ++child)
389  if (cell->child(child)->get_fe().dofs_per_cell >
390  cell->child(most_general_child)->get_fe().dofs_per_cell)
391  most_general_child = child;
392  const unsigned int target_fe_index =
393  cell->child(most_general_child)->active_fe_index();
394 
395  const unsigned int dofs_per_cell =
396  dof_handler->get_fe(target_fe_index).dofs_per_cell;
397 
398  std::vector<Vector<typename VectorType::value_type>>(
399  in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
400  .swap(dof_values_on_cell[n_cf]);
401 
402 
403  // store the data of each of the input vectors. get this data
404  // as interpolated onto a finite element space that encompasses
405  // that of all the children. note that
406  // cell->get_interpolated_dof_values already does all of the
407  // interpolations between spaces
408  for (unsigned int j = 0; j < in_size; ++j)
409  cell->get_interpolated_dof_values(all_in[j],
410  dof_values_on_cell[n_cf][j],
411  target_fe_index);
412  cell_map[std::make_pair(cell->level(), cell->index())] =
413  Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
414  ++n_cf;
415  }
416  }
417  Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
418  Assert(n_cf == n_coarsen_fathers, ExcInternalError());
419 
421 }
422 
423 
424 
425 template <int dim, typename VectorType, typename DoFHandlerType>
426 void
429 {
430  std::vector<VectorType> all_in(1, in);
432 }
433 
434 
435 
436 template <int dim, typename VectorType, typename DoFHandlerType>
437 void
439  const std::vector<VectorType> &all_in,
440  std::vector<VectorType> & all_out) const
441 {
443  const unsigned int size = all_in.size();
444  Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
445  for (unsigned int i = 0; i < size; ++i)
446  Assert(all_in[i].size() == n_dofs_old,
447  ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
448  for (unsigned int i = 0; i < all_out.size(); ++i)
449  Assert(all_out[i].size() == dof_handler->n_dofs(),
450  ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
451  for (unsigned int i = 0; i < size; ++i)
452  for (unsigned int j = 0; j < size; ++j)
453  Assert(&all_in[i] != &all_out[j],
454  ExcMessage("Vectors cannot be used as input and output"
455  " at the same time!"));
456 
458  std::vector<types::global_dof_index> dofs;
459 
460  typename std::map<std::pair<unsigned int, unsigned int>,
461  Pointerstruct>::const_iterator pointerstruct,
462  cell_map_end = cell_map.end();
463 
464  Table<2, FullMatrix<double>> interpolation_hp;
467 
468  typename DoFHandlerType::cell_iterator cell = dof_handler->begin(),
469  endc = dof_handler->end();
470  for (; cell != endc; ++cell)
471  {
472  pointerstruct =
473  cell_map.find(std::make_pair(cell->level(), cell->index()));
474 
475  if (pointerstruct != cell_map_end)
476  {
477  const std::vector<types::global_dof_index> *const indexptr =
478  pointerstruct->second.indices_ptr;
479 
480  const std::vector<Vector<typename VectorType::value_type>>
481  *const valuesptr = pointerstruct->second.dof_values_ptr;
482 
483  // cell stayed as it was or was refined
484  if (indexptr)
485  {
486  Assert(valuesptr == nullptr, ExcInternalError());
487 
488  const unsigned int old_fe_index =
489  pointerstruct->second.active_fe_index;
490 
491  // get the values of
492  // each of the input
493  // data vectors on this
494  // cell and prolong it
495  // to its children
496  unsigned int in_size = indexptr->size();
497  for (unsigned int j = 0; j < size; ++j)
498  {
499  tmp.reinit(in_size, true);
500  for (unsigned int i = 0; i < in_size; ++i)
501  tmp(i) =
502  internal::ElementAccess<VectorType>::get(all_in[j],
503  (*indexptr)[i]);
504 
505  cell->set_dof_values_by_interpolation(tmp,
506  all_out[j],
507  old_fe_index);
508  }
509  }
510  else if (valuesptr)
511  // the children of this cell were
512  // deleted
513  {
514  Assert(!cell->has_children(), ExcInternalError());
515  Assert(indexptr == nullptr, ExcInternalError());
516 
517  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
518  dofs.resize(dofs_per_cell);
519  // get the local
520  // indices
521  cell->get_dof_indices(dofs);
522 
523  // distribute the
524  // stored data to the
525  // new vectors
526  for (unsigned int j = 0; j < size; ++j)
527  {
528  // make sure that the size of
529  // the stored indices is the
530  // same as
531  // dofs_per_cell. this is
532  // kind of a test if we use
533  // the same fe in the hp
534  // case. to really do that
535  // test we would have to
536  // store the fe_index of all
537  // cells
538  const Vector<typename VectorType::value_type> *data = nullptr;
539  const unsigned int active_fe_index = cell->active_fe_index();
540  if (active_fe_index != pointerstruct->second.active_fe_index)
541  {
542  const unsigned int old_index =
543  pointerstruct->second.active_fe_index;
544  tmp.reinit(dofs_per_cell, true);
546  (*valuesptr)[j].size(),
547  interpolation_hp(active_fe_index, old_index).n());
549  tmp.size(),
550  interpolation_hp(active_fe_index, old_index).m());
551  interpolation_hp(active_fe_index, old_index)
552  .vmult(tmp, (*valuesptr)[j]);
553  data = &tmp;
554  }
555  else
556  data = &(*valuesptr)[j];
557 
558 
559  for (unsigned int i = 0; i < dofs_per_cell; ++i)
561  dofs[i],
562  all_out[j]);
563  }
564  }
565  // undefined status
566  else
567  Assert(false, ExcInternalError());
568  }
569  }
570 }
571 
572 
573 
574 template <int dim, typename VectorType, typename DoFHandlerType>
575 void
577  const VectorType &in,
578  VectorType & out) const
579 {
580  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
581  Assert(out.size() == dof_handler->n_dofs(),
582  ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
583 
584  std::vector<VectorType> all_in(1);
585  all_in[0] = in;
586  std::vector<VectorType> all_out(1);
587  all_out[0] = out;
588  interpolate(all_in, all_out);
589  out = all_out[0];
590 }
591 
592 
593 
594 template <int dim, typename VectorType, typename DoFHandlerType>
595 std::size_t
597 {
598  // at the moment we do not include the memory
599  // consumption of the cell_map as we have no
600  // real idea about memory consumption of a
601  // std::map
604  sizeof(prepared_for) +
607 }
608 
609 
610 
611 template <int dim, typename VectorType, typename DoFHandlerType>
612 std::size_t
614  memory_consumption() const
615 {
616  return sizeof(*this);
617 }
618 
619 
620 /*-------------- Explicit Instantiations -------------------------------*/
621 #define SPLIT_INSTANTIATIONS_COUNT 4
622 #ifndef SPLIT_INSTANTIATIONS_INDEX
623 # define SPLIT_INSTANTIATIONS_INDEX 0
624 #endif
625 #include "solution_transfer.inst"
626 
627 DEAL_II_NAMESPACE_CLOSE
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
types::global_dof_index n_dofs_old
SmartPointer< const DoFHandlerType, SolutionTransfer< dim, VectorType, DoFHandlerType > > dof_handler
static ::ExceptionBase & ExcAlreadyPrepForRef()
static ::ExceptionBase & ExcNotPrepared()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::map< std::pair< unsigned int, unsigned int >, Pointerstruct > cell_map
std::vector< std::vector< types::global_dof_index > > indices_on_cell
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void prepare_for_pure_refinement()
static ::ExceptionBase & ExcAlreadyPrepForCoarseAndRef()
void swap(SmartPointer< T, Q > &tt)
Definition: smartpointer.h:391
void extract_interpolation_matrices(const DoFHandlerType &, ::Table< 2, FullMatrix< double >> &)
SolutionTransfer(const DoFHandlerType &dof)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
__global__ void set(Number *val, const Number s, const size_type N)
void prepare_for_coarsening_and_refinement(const std::vector< VectorType > &all_in)
void refine_interpolate(const VectorType &in, VectorType &out) const
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
Definition: table.h:37
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)
static ::ExceptionBase & ExcInternalError()