Reference documentation for deal.II version Git 713825e468 2021-05-17 16:05:53 -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) 2009 - 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 
16 
17 #include <deal.II/base/config.h>
18 
19 #ifdef DEAL_II_WITH_P4EST
20 
23 
25 # include <deal.II/dofs/dof_tools.h>
26 
29 
30 # include <deal.II/hp/dof_handler.h>
31 
39 # include <deal.II/lac/vector.h>
40 
41 # include <functional>
42 # include <numeric>
43 
44 
46 
47 
48 namespace
49 {
58  template <typename value_type>
59  std::vector<char>
60  pack_dof_values(std::vector<Vector<value_type>> &dof_values,
61  const unsigned int dofs_per_cell)
62  {
63  for (const auto &values : dof_values)
64  {
65  AssertDimension(values.size(), dofs_per_cell);
66  (void)values;
67  }
68 
69  const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
70 
71  std::vector<char> buffer(dof_values.size() * bytes_per_entry);
72  for (unsigned int i = 0; i < dof_values.size(); ++i)
73  std::memcpy(&buffer[i * bytes_per_entry],
74  &dof_values[i](0),
75  bytes_per_entry);
76 
77  return buffer;
78  }
79 
80 
81 
85  template <typename value_type>
86  std::vector<Vector<value_type>>
87  unpack_dof_values(
88  const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
89  const unsigned int dofs_per_cell)
90  {
91  const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
92  const unsigned int n_elements = data_range.size() / bytes_per_entry;
93 
94  Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
95 
96  std::vector<Vector<value_type>> unpacked_data;
97  unpacked_data.reserve(n_elements);
98  for (unsigned int i = 0; i < n_elements; ++i)
99  {
100  Vector<value_type> dof_values(dofs_per_cell);
101  std::memcpy(&dof_values(0),
102  &(*std::next(data_range.begin(), i * bytes_per_entry)),
103  bytes_per_entry);
104  unpacked_data.emplace_back(std::move(dof_values));
105  }
106 
107  return unpacked_data;
108  }
109 } // namespace
110 
111 
112 
113 namespace parallel
114 {
115  namespace distributed
116  {
117  template <int dim, typename VectorType, typename DoFHandlerType>
119  const DoFHandlerType &dof)
120  : dof_handler(&dof, typeid(*this).name())
121  , handle(numbers::invalid_unsigned_int)
122  {
123  Assert(
124  (dynamic_cast<const parallel::DistributedTriangulationBase<
125  dim,
126  DoFHandlerType::space_dimension> *>(
127  &dof_handler->get_triangulation()) != nullptr),
128  ExcMessage(
129  "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
130  }
131 
132 
133 
134  template <int dim, typename VectorType, typename DoFHandlerType>
135  void
138  const std::vector<const VectorType *> &all_in)
139  {
140  for (unsigned int i = 0; i < all_in.size(); ++i)
141  Assert(all_in[i]->size() == dof_handler->n_dofs(),
142  ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
143 
144  input_vectors = all_in;
146  }
147 
148 
149 
150  template <int dim, typename VectorType, typename DoFHandlerType>
151  void
153  {
154  // TODO: casting away constness is bad
155  auto *tria = (dynamic_cast<parallel::DistributedTriangulationBase<
156  dim,
157  DoFHandlerType::space_dimension> *>(
159  *>(&dof_handler->get_triangulation())));
160  Assert(tria != nullptr, ExcInternalError());
161 
162  handle = tria->register_data_attach(
163  [this](
165  cell_iterator &cell_,
167  CellStatus status) { return this->pack_callback(cell_, status); },
168  /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
169  }
170 
171 
172 
173  template <int dim, typename VectorType, typename DoFHandlerType>
174  void
177  {
178  std::vector<const VectorType *> all_in(1, &in);
180  }
181 
182 
183 
184  template <int dim, typename VectorType, typename DoFHandlerType>
185  void
188  {
189  std::vector<const VectorType *> all_in(1, &in);
191  }
192 
193 
194 
195  template <int dim, typename VectorType, typename DoFHandlerType>
196  void
198  prepare_for_serialization(const std::vector<const VectorType *> &all_in)
199  {
201  }
202 
203 
204 
205  template <int dim, typename VectorType, typename DoFHandlerType>
206  void
208  VectorType &in)
209  {
210  std::vector<VectorType *> all_in(1, &in);
211  deserialize(all_in);
212  }
213 
214 
215 
216  template <int dim, typename VectorType, typename DoFHandlerType>
217  void
219  std::vector<VectorType *> &all_in)
220  {
222 
223  // this makes interpolate() happy
224  input_vectors.resize(all_in.size());
225 
226  interpolate(all_in);
227  }
228 
229 
230  template <int dim, typename VectorType, typename DoFHandlerType>
231  void
233  std::vector<VectorType *> &all_out)
234  {
235  Assert(input_vectors.size() == all_out.size(),
236  ExcDimensionMismatch(input_vectors.size(), all_out.size()));
237  for (unsigned int i = 0; i < all_out.size(); ++i)
238  Assert(all_out[i]->size() == dof_handler->n_dofs(),
239  ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
240 
241  // TODO: casting away constness is bad
242  auto *tria = (dynamic_cast<parallel::DistributedTriangulationBase<
243  dim,
244  DoFHandlerType::space_dimension> *>(
246  *>(&dof_handler->get_triangulation())));
247  Assert(tria != nullptr, ExcInternalError());
248 
249  tria->notify_ready_to_unpack(
250  handle,
251  [this, &all_out](
253  cell_iterator &cell_,
255  CellStatus status,
256  const boost::iterator_range<std::vector<char>::const_iterator>
257  &data_range) {
258  this->unpack_callback(cell_, status, data_range, all_out);
259  });
260 
261  for (typename std::vector<VectorType *>::iterator it = all_out.begin();
262  it != all_out.end();
263  ++it)
264  (*it)->compress(::VectorOperation::insert);
265 
266  input_vectors.clear();
267  }
268 
269 
270 
271  template <int dim, typename VectorType, typename DoFHandlerType>
272  void
274  VectorType &out)
275  {
276  std::vector<VectorType *> all_out(1, &out);
277  interpolate(all_out);
278  }
279 
280 
281 
282  template <int dim, typename VectorType, typename DoFHandlerType>
283  std::vector<char>
286  cell_iterator &cell_,
287  const typename Triangulation<dim,
288  DoFHandlerType::space_dimension>::CellStatus
289  status)
290  {
291  typename DoFHandlerType::cell_iterator cell(*cell_, dof_handler);
292 
293  // create buffer for each individual object
294  std::vector<::Vector<typename VectorType::value_type>> dof_values(
295  input_vectors.size());
296 
297  unsigned int fe_index = 0;
298  if (dof_handler->has_hp_capabilities())
299  {
300  switch (status)
301  {
303  dim,
304  DoFHandlerType::space_dimension>::CELL_PERSIST:
306  dim,
307  DoFHandlerType::space_dimension>::CELL_REFINE:
308  {
309  fe_index = cell->future_fe_index();
310  break;
311  }
312 
314  dim,
315  DoFHandlerType::space_dimension>::CELL_COARSEN:
316  {
317  // In case of coarsening, we need to find a suitable FE index
318  // for the parent cell. We choose the 'least dominant fe'
319  // on all children from the associated FECollection.
320 # ifdef DEBUG
321  for (const auto &child : cell->child_iterators())
322  Assert(child->is_active() && child->coarsen_flag_set(),
324  dim>::ExcInconsistentCoarseningFlags());
325 # endif
326 
327  fe_index = ::internal::hp::DoFHandlerImplementation::
328  dominated_future_fe_on_children<
329  dim,
330  DoFHandlerType::space_dimension>(cell);
331  break;
332  }
333 
334  default:
335  Assert(false, ExcInternalError());
336  break;
337  }
338  }
339 
340  const unsigned int dofs_per_cell =
341  dof_handler->get_fe(fe_index).n_dofs_per_cell();
342 
343  if (dofs_per_cell == 0)
344  return std::vector<char>(); // nothing to do for FE_Nothing
345 
346  auto it_input = input_vectors.cbegin();
347  auto it_output = dof_values.begin();
348  for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
349  {
350  it_output->reinit(dofs_per_cell);
351  cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
352  }
353 
354  return pack_dof_values<typename VectorType::value_type>(dof_values,
355  dofs_per_cell);
356  }
357 
358 
359 
360  template <int dim, typename VectorType, typename DoFHandlerType>
361  void
364  cell_iterator &cell_,
365  const typename Triangulation<dim,
366  DoFHandlerType::space_dimension>::CellStatus
367  status,
368  const boost::iterator_range<std::vector<char>::const_iterator>
369  & data_range,
370  std::vector<VectorType *> &all_out)
371  {
372  typename DoFHandlerType::cell_iterator cell(*cell_, dof_handler);
373 
374  unsigned int fe_index = 0;
375  if (dof_handler->has_hp_capabilities())
376  {
377  switch (status)
378  {
380  dim,
381  DoFHandlerType::space_dimension>::CELL_PERSIST:
383  dim,
384  DoFHandlerType::space_dimension>::CELL_COARSEN:
385  {
386  fe_index = cell->active_fe_index();
387  break;
388  }
389 
391  dim,
392  DoFHandlerType::space_dimension>::CELL_REFINE:
393  {
394  // After refinement, this particular cell is no longer active,
395  // and its children have inherited its FE index. However, to
396  // unpack the data on the old cell, we need to recover its FE
397  // index from one of the children. Just to be sure, we also
398  // check if all children have the same FE index.
399  fe_index = cell->child(0)->active_fe_index();
400  for (unsigned int child_index = 1;
401  child_index < cell->n_children();
402  ++child_index)
403  Assert(cell->child(child_index)->active_fe_index() ==
404  fe_index,
405  ExcInternalError());
406  break;
407  }
408 
409  default:
410  Assert(false, ExcInternalError());
411  break;
412  }
413  }
414 
415  const unsigned int dofs_per_cell =
416  dof_handler->get_fe(fe_index).n_dofs_per_cell();
417 
418  if (dofs_per_cell == 0)
419  return; // nothing to do for FE_Nothing
420 
421  const std::vector<::Vector<typename VectorType::value_type>>
422  dof_values =
423  unpack_dof_values<typename VectorType::value_type>(data_range,
424  dofs_per_cell);
425 
426  // check if sizes match
427  Assert(dof_values.size() == all_out.size(), ExcInternalError());
428 
429  // check if we have enough dofs provided by the FE object
430  // to interpolate the transferred data correctly
431  for (auto it_dof_values = dof_values.begin();
432  it_dof_values != dof_values.end();
433  ++it_dof_values)
434  Assert(
435  dofs_per_cell == it_dof_values->size(),
436  ExcMessage(
437  "The transferred data was packed with a different number of dofs than the "
438  "currently registered FE object assigned to the DoFHandler has."));
439 
440  // distribute data for each registered vector on mesh
441  auto it_input = dof_values.cbegin();
442  auto it_output = all_out.begin();
443  for (; it_input != dof_values.cend(); ++it_input, ++it_output)
444  cell->set_dof_values_by_interpolation(*it_input,
445  *(*it_output),
446  fe_index);
447  }
448  } // namespace distributed
449 } // namespace parallel
450 
451 
452 // explicit instantiations
453 # include "solution_transfer.inst"
454 
456 
457 #endif
void prepare_for_coarsening_and_refinement(const std::vector< const VectorType *> &all_in)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
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:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:396
std::vector< const VectorType * > input_vectors
SmartPointer< const DoFHandlerType, SolutionTransfer< dim, VectorType, DoFHandlerType > > dof_handler
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
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()