Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solution_transfer.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 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 
17 #include <deal.II/base/config.h>
18 
19 #ifdef DEAL_II_WITH_P4EST
20 
21 # include <deal.II/distributed/solution_transfer.h>
22 # include <deal.II/distributed/tria.h>
23 
24 # include <deal.II/dofs/dof_accessor.h>
25 # include <deal.II/dofs/dof_tools.h>
26 
27 # include <deal.II/grid/tria_accessor.h>
28 # include <deal.II/grid/tria_iterator.h>
29 
30 # include <deal.II/hp/dof_handler.h>
31 
32 # include <deal.II/lac/block_vector.h>
33 # include <deal.II/lac/la_parallel_block_vector.h>
34 # include <deal.II/lac/la_parallel_vector.h>
35 # include <deal.II/lac/petsc_block_vector.h>
36 # include <deal.II/lac/petsc_vector.h>
37 # include <deal.II/lac/trilinos_parallel_block_vector.h>
38 # include <deal.II/lac/trilinos_vector.h>
39 # include <deal.II/lac/vector.h>
40 
41 # include <functional>
42 
43 DEAL_II_NAMESPACE_OPEN
44 
45 namespace parallel
46 {
47  namespace distributed
48  {
49  template <int dim, typename VectorType, typename DoFHandlerType>
51  const DoFHandlerType &dof)
52  : dof_handler(&dof, typeid(*this).name())
53  , handle(numbers::invalid_unsigned_int)
54  {
55  Assert(
56  (dynamic_cast<const parallel::distributed::
57  Triangulation<dim, DoFHandlerType::space_dimension> *>(
58  &dof_handler->get_triangulation()) != nullptr),
59  ExcMessage(
60  "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
61  }
62 
63 
64 
65  template <int dim, typename VectorType, typename DoFHandlerType>
66  void
69  const std::vector<const VectorType *> &all_in)
70  {
71  input_vectors = all_in;
72  register_data_attach();
73  }
74 
75 
76 
77  template <int dim, typename VectorType, typename DoFHandlerType>
78  void
80  {
81  // TODO: casting away constness is bad
83  *tria = (dynamic_cast<parallel::distributed::Triangulation<
84  dim,
85  DoFHandlerType::space_dimension> *>(
87  *>(&dof_handler->get_triangulation())));
88  Assert(tria != nullptr, ExcInternalError());
89 
90  handle = tria->register_data_attach(
91  std::bind(
93  this,
94  std::placeholders::_1,
95  std::placeholders::_2),
96  /*returns_variable_size_data=*/DoFHandlerType::is_hp_dof_handler);
97  }
98 
99 
100 
101  template <int dim, typename VectorType, typename DoFHandlerType>
102  void
105  {
106  std::vector<const VectorType *> all_in(1, &in);
107  prepare_for_coarsening_and_refinement(all_in);
108  }
109 
110 
111 
112  template <int dim, typename VectorType, typename DoFHandlerType>
113  void
115  prepare_for_serialization(const VectorType &in)
116  {
117  std::vector<const VectorType *> all_in(1, &in);
118  prepare_for_serialization(all_in);
119  }
120 
121 
122 
123  template <int dim, typename VectorType, typename DoFHandlerType>
124  void
126  prepare_for_serialization(const std::vector<const VectorType *> &all_in)
127  {
128  prepare_for_coarsening_and_refinement(all_in);
129  }
130 
131 
132 
133  template <int dim, typename VectorType, typename DoFHandlerType>
134  void
136  const VectorType &in)
137  {
138  prepare_for_serialization(in);
139  }
140 
141 
142 
143  template <int dim, typename VectorType, typename DoFHandlerType>
144  void
146  const std::vector<const VectorType *> &all_in)
147  {
148  prepare_for_serialization(all_in);
149  }
150 
151 
152 
153  template <int dim, typename VectorType, typename DoFHandlerType>
154  void
156  VectorType &in)
157  {
158  std::vector<VectorType *> all_in(1, &in);
159  deserialize(all_in);
160  }
161 
162 
163 
164  template <int dim, typename VectorType, typename DoFHandlerType>
165  void
167  std::vector<VectorType *> &all_in)
168  {
169  register_data_attach();
170 
171  // this makes interpolate() happy
172  input_vectors.resize(all_in.size());
173 
174  interpolate(all_in);
175  }
176 
177 
178  template <int dim, typename VectorType, typename DoFHandlerType>
179  void
181  std::vector<VectorType *> &all_out)
182  {
183  Assert(input_vectors.size() == all_out.size(),
184  ExcDimensionMismatch(input_vectors.size(), all_out.size()));
185 
186  // TODO: casting away constness is bad
188  *tria = (dynamic_cast<parallel::distributed::Triangulation<
189  dim,
190  DoFHandlerType::space_dimension> *>(
192  *>(&dof_handler->get_triangulation())));
193  Assert(tria != nullptr, ExcInternalError());
194 
196  handle,
197  std::bind(
199  this,
200  std::placeholders::_1,
201  std::placeholders::_2,
202  std::placeholders::_3,
203  std::ref(all_out)));
204 
205 
206  for (typename std::vector<VectorType *>::iterator it = all_out.begin();
207  it != all_out.end();
208  ++it)
209  (*it)->compress(::VectorOperation::insert);
210 
211  input_vectors.clear();
212  }
213 
214 
215 
216  template <int dim, typename VectorType, typename DoFHandlerType>
217  void
219  VectorType &out)
220  {
221  std::vector<VectorType *> all_out(1, &out);
222  interpolate(all_out);
223  }
224 
225 
226 
227  template <int dim, typename VectorType, typename DoFHandlerType>
228  std::vector<char>
231  cell_iterator &cell_,
232  const typename Triangulation<dim,
233  DoFHandlerType::space_dimension>::CellStatus
234  status)
235  {
236  typename DoFHandlerType::cell_iterator cell(*cell_, dof_handler);
237 
238  // create buffer for each individual object
239  std::vector<::Vector<typename VectorType::value_type>> dofvalues(
240  input_vectors.size());
241 
242  if (DoFHandlerType::is_hp_dof_handler)
243  {
244  unsigned int fe_index = numbers::invalid_unsigned_int;
245  switch (status)
246  {
248  dim,
249  DoFHandlerType::space_dimension>::CELL_PERSIST:
251  dim,
252  DoFHandlerType::space_dimension>::CELL_REFINE:
253  {
254  fe_index = cell->future_fe_index();
255  break;
256  }
257 
259  dim,
260  DoFHandlerType::space_dimension>::CELL_COARSEN:
261  {
262  // In case of coarsening, we need to find a suitable fe index
263  // for the parent cell. We choose the 'least dominating fe'
264  // on all children from the associated FECollection.
265  std::set<unsigned int> fe_indices_children;
266  for (unsigned int child_index = 0;
267  child_index < GeometryInfo<dim>::max_children_per_cell;
268  ++child_index)
269  fe_indices_children.insert(
270  cell->child(child_index)->future_fe_index());
271 
272  fe_index = dof_handler->get_fe_collection()
273  .find_dominating_fe_extended(fe_indices_children,
274  /*codim=*/0);
275 
276  Assert(
277  fe_index != numbers::invalid_unsigned_int,
278  ExcMessage(
279  "No FiniteElement has been found in your FECollection "
280  "that dominates all children of a cell you are trying "
281  "to coarsen!"));
282 
283  break;
284  }
285 
286  default:
287  Assert(false, ExcInternalError());
288  break;
289  }
290 
291  const unsigned int dofs_per_cell =
292  dof_handler->get_fe(fe_index).dofs_per_cell;
293 
294  auto it_input = input_vectors.cbegin();
295  auto it_output = dofvalues.begin();
296  for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
297  {
298  it_output->reinit(dofs_per_cell);
299  cell->get_interpolated_dof_values(*(*it_input),
300  *it_output,
301  fe_index);
302  }
303  }
304  else
305  {
306  auto it_input = input_vectors.cbegin();
307  auto it_output = dofvalues.begin();
308  for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
309  {
310  it_output->reinit(cell->get_fe().dofs_per_cell);
311  cell->get_interpolated_dof_values(*(*it_input), *it_output);
312  }
313  }
314 
315  // We choose to compress the packed data depending on the registered
316  // dofhandler object. In case of hp, we write in the variable size buffer
317  // and thus allow compression. In the other case, we use the fixed size
318  // buffer and require consistent data sizes on each cell, so we leave it
319  // uncompressed.
320  return Utilities::pack(
321  dofvalues, /*allow_compression=*/DoFHandlerType::is_hp_dof_handler);
322  }
323 
324 
325 
326  template <int dim, typename VectorType, typename DoFHandlerType>
327  void
330  cell_iterator &cell_,
331  const typename Triangulation<dim,
332  DoFHandlerType::space_dimension>::CellStatus
333  status,
334  const boost::iterator_range<std::vector<char>::const_iterator>
335  & data_range,
336  std::vector<VectorType *> &all_out)
337  {
338  typename DoFHandlerType::cell_iterator cell(*cell_, dof_handler);
339 
340  const std::vector<::Vector<typename VectorType::value_type>>
341  dofvalues = Utilities::unpack<
342  std::vector<::Vector<typename VectorType::value_type>>>(
343  data_range.begin(),
344  data_range.end(),
345  /*allow_compression=*/DoFHandlerType::is_hp_dof_handler);
346 
347  // check if sizes match
348  Assert(dofvalues.size() == all_out.size(), ExcInternalError());
349 
350  if (DoFHandlerType::is_hp_dof_handler)
351  {
352  unsigned int fe_index = numbers::invalid_unsigned_int;
353 
354  switch (status)
355  {
357  dim,
358  DoFHandlerType::space_dimension>::CELL_PERSIST:
360  dim,
361  DoFHandlerType::space_dimension>::CELL_COARSEN:
362  {
363  fe_index = cell->active_fe_index();
364  break;
365  }
366 
368  dim,
369  DoFHandlerType::space_dimension>::CELL_REFINE:
370  {
371  // After refinement, this particular cell is no longer active,
372  // and its children have inherited its fe index. However, to
373  // unpack the data on the old cell, we need to recover its fe
374  // index from one of the children. Just to be sure, we also
375  // check if all children have the same fe index.
376  fe_index = cell->child(0)->active_fe_index();
377  for (unsigned int child_index = 1;
378  child_index < GeometryInfo<dim>::max_children_per_cell;
379  ++child_index)
380  Assert(cell->child(child_index)->active_fe_index() ==
381  fe_index,
382  ExcInternalError());
383  break;
384  }
385 
386  default:
387  Assert(false, ExcInternalError());
388  break;
389  }
390 
391  // check if we have enough dofs provided by the FE object
392  // to interpolate the transferred data correctly
393  for (auto it_dofvalues = dofvalues.begin();
394  it_dofvalues != dofvalues.end();
395  ++it_dofvalues)
396  Assert(
397  dof_handler->get_fe(fe_index).dofs_per_cell ==
398  it_dofvalues->size(),
399  ExcMessage(
400  "The transferred data was packed with a different number of dofs than the"
401  "currently registered FE object assigned to the DoFHandler has."));
402 
403  // distribute data for each registered vector on mesh
404  auto it_input = dofvalues.cbegin();
405  auto it_output = all_out.begin();
406  for (; it_input != dofvalues.cend(); ++it_input, ++it_output)
407  cell->set_dof_values_by_interpolation(*it_input,
408  *(*it_output),
409  fe_index);
410  }
411  else
412  {
413  // check if we have enough dofs provided by the FE object
414  // to interpolate the transferred data correctly
415  for (auto it_dofvalues = dofvalues.begin();
416  it_dofvalues != dofvalues.end();
417  ++it_dofvalues)
418  Assert(
419  cell->get_fe().dofs_per_cell == it_dofvalues->size(),
420  ExcMessage(
421  "The transferred data was packed with a different number of dofs than the"
422  "currently registered FE object assigned to the DoFHandler has."));
423 
424  // distribute data for each registered vector on mesh
425  auto it_input = dofvalues.cbegin();
426  auto it_output = all_out.begin();
427  for (; it_input != dofvalues.cend(); ++it_input, ++it_output)
428  cell->set_dof_values_by_interpolation(*it_input, *(*it_output));
429  }
430  }
431 
432 
433  } // namespace distributed
434 } // namespace parallel
435 
436 
437 // explicit instantiations
438 # include "solution_transfer.inst"
439 
440 DEAL_II_NAMESPACE_CLOSE
441 
442 #endif
void prepare_for_coarsening_and_refinement(const std::vector< const VectorType *> &all_in)
unsigned int register_data_attach(const std::function< std::vector< char >(const cell_iterator &, const CellStatus)> &pack_callback, const bool returns_variable_size_data)
Definition: tria.cc:4357
static const unsigned int invalid_unsigned_int
Definition: types.h:173
std::vector< char > pack_callback(const typename Triangulation< dim, DoFHandlerType::space_dimension >::cell_iterator &cell, const typename Triangulation< dim, DoFHandlerType::space_dimension >::CellStatus status)
SolutionTransfer(const DoFHandlerType &dof)
void interpolate(std::vector< VectorType *> &all_out)
void prepare_for_serialization(const VectorType &in)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void notify_ready_to_unpack(const unsigned int handle, const std::function< void(const cell_iterator &, const CellStatus, const boost::iterator_range< std::vector< char >::const_iterator > &)> &unpack_callback)
Definition: tria.cc:4387
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1174
void prepare_serialization(const VectorType &in)
SmartPointer< const DoFHandlerType, SolutionTransfer< dim, VectorType, DoFHandlerType > > dof_handler
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1318
void unpack_callback(const typename Triangulation< dim, DoFHandlerType::space_dimension >::cell_iterator &cell, const typename Triangulation< dim, DoFHandlerType::space_dimension >::CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range, std::vector< VectorType *> &all_out)
static ::ExceptionBase & ExcInternalError()