Reference documentation for deal.II version GIT c545eda070 2023-01-27 00:25:02+00:00
\(\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 - 2022 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 
37 # include <deal.II/lac/vector.h>
38 
39 # include <functional>
40 # include <numeric>
41 
42 
44 
45 
46 namespace
47 {
56  template <typename value_type>
57  std::vector<char>
58  pack_dof_values(std::vector<Vector<value_type>> &dof_values,
59  const unsigned int dofs_per_cell)
60  {
61  for (const auto &values : dof_values)
62  {
63  AssertDimension(values.size(), dofs_per_cell);
64  (void)values;
65  }
66 
67  const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
68 
69  std::vector<char> buffer(dof_values.size() * bytes_per_entry);
70  for (unsigned int i = 0; i < dof_values.size(); ++i)
71  std::memcpy(&buffer[i * bytes_per_entry],
72  &dof_values[i](0),
73  bytes_per_entry);
74 
75  return buffer;
76  }
77 
78 
79 
83  template <typename value_type>
84  std::vector<Vector<value_type>>
85  unpack_dof_values(
86  const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
87  const unsigned int dofs_per_cell)
88  {
89  const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
90  const unsigned int n_elements = data_range.size() / bytes_per_entry;
91 
92  Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
93 
94  std::vector<Vector<value_type>> unpacked_data;
95  unpacked_data.reserve(n_elements);
96  for (unsigned int i = 0; i < n_elements; ++i)
97  {
98  Vector<value_type> dof_values(dofs_per_cell);
99  std::memcpy(&dof_values(0),
100  &(*std::next(data_range.begin(), i * bytes_per_entry)),
101  bytes_per_entry);
102  unpacked_data.emplace_back(std::move(dof_values));
103  }
104 
105  return unpacked_data;
106  }
107 } // namespace
108 
109 
110 
111 namespace parallel
112 {
113  namespace distributed
114  {
115  template <int dim, typename VectorType, int spacedim>
117  const DoFHandler<dim, spacedim> &dof,
118  const bool average_values)
119  : dof_handler(&dof, typeid(*this).name())
120  , average_values(average_values)
121  , handle(numbers::invalid_unsigned_int)
122  {
123  Assert(
124  (dynamic_cast<
126  &dof_handler->get_triangulation()) != nullptr),
127  ExcMessage(
128  "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
129  }
130 
131 
132 
133  template <int dim, typename VectorType, int spacedim>
134  void
137  const std::vector<const VectorType *> &all_in)
138  {
139  for (unsigned int i = 0; i < all_in.size(); ++i)
140  Assert(all_in[i]->size() == dof_handler->n_dofs(),
141  ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
142 
143  input_vectors = all_in;
144  register_data_attach();
145  }
146 
147 
148 
149  template <int dim, typename VectorType, int spacedim>
150  void
152  {
153  // TODO: casting away constness is bad
156  const_cast<::Triangulation<dim, spacedim> *>(
157  &dof_handler->get_triangulation())));
158  Assert(tria != nullptr, ExcInternalError());
159 
161  ExcMessage("You can only add one solution per "
162  "SolutionTransfer object."));
163 
164  handle = tria->register_data_attach(
165  [this](
166  const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
167  const typename Triangulation<dim, spacedim>::CellStatus status) {
168  return this->pack_callback(cell_, status);
169  },
170  /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
171  }
172 
173 
174 
175  template <int dim, typename VectorType, int spacedim>
176  void
178  prepare_for_coarsening_and_refinement(const VectorType &in)
179  {
180  std::vector<const VectorType *> all_in(1, &in);
181  prepare_for_coarsening_and_refinement(all_in);
182  }
183 
184 
185 
186  template <int dim, typename VectorType, int spacedim>
187  void
189  const VectorType &in)
190  {
191  std::vector<const VectorType *> all_in(1, &in);
192  prepare_for_serialization(all_in);
193  }
194 
195 
196 
197  template <int dim, typename VectorType, int spacedim>
198  void
200  const std::vector<const VectorType *> &all_in)
201  {
202  prepare_for_coarsening_and_refinement(all_in);
203  }
204 
205 
206 
207  template <int dim, typename VectorType, int spacedim>
208  void
210  {
211  std::vector<VectorType *> all_in(1, &in);
212  deserialize(all_in);
213  }
214 
215 
216 
217  template <int dim, typename VectorType, int spacedim>
218  void
220  std::vector<VectorType *> &all_in)
221  {
222  register_data_attach();
223 
224  // this makes interpolate() happy
225  input_vectors.resize(all_in.size());
226 
227  interpolate(all_in);
228  }
229 
230 
231  template <int dim, typename VectorType, int spacedim>
232  void
234  std::vector<VectorType *> &all_out)
235  {
236  Assert(input_vectors.size() == all_out.size(),
237  ExcDimensionMismatch(input_vectors.size(), all_out.size()));
238  for (unsigned int i = 0; i < all_out.size(); ++i)
239  Assert(all_out[i]->size() == dof_handler->n_dofs(),
240  ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
241 
242  // TODO: casting away constness is bad
245  const_cast<::Triangulation<dim, spacedim> *>(
246  &dof_handler->get_triangulation())));
247  Assert(tria != nullptr, ExcInternalError());
248 
249  if (average_values)
250  for (const auto vec : all_out)
251  *vec = 0.0;
252 
253  VectorType valence;
254 
255  // initialize valence vector only if we need to average
256  if (average_values)
257  valence.reinit(*all_out[0]);
258 
259  tria->notify_ready_to_unpack(
260  handle,
261  [this, &all_out, &valence](
262  const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
263  const typename Triangulation<dim, spacedim>::CellStatus status,
264  const boost::iterator_range<std::vector<char>::const_iterator>
265  &data_range) {
266  this->unpack_callback(cell_, status, data_range, all_out, valence);
267  });
268 
269  if (average_values)
270  {
271  // finalize valence: compress and invert
272  using Number = typename VectorType::value_type;
273  valence.compress(::VectorOperation::add);
274  for (const auto i : valence.locally_owned_elements())
275  valence[i] = (static_cast<Number>(valence[i]) == Number() ?
276  Number() :
277  (Number(1.0) / static_cast<Number>(valence[i])));
278  valence.compress(::VectorOperation::insert);
279 
280  for (const auto vec : all_out)
281  {
282  // compress and weight with valence
283  vec->compress(::VectorOperation::add);
284  vec->scale(valence);
285  }
286  }
287  else
288  {
289  for (const auto vec : all_out)
290  vec->compress(::VectorOperation::insert);
291  }
292 
293  input_vectors.clear();
295  }
296 
297 
298 
299  template <int dim, typename VectorType, int spacedim>
300  void
302  {
303  std::vector<VectorType *> all_out(1, &out);
304  interpolate(all_out);
305  }
306 
307 
308 
309  template <int dim, typename VectorType, int spacedim>
310  std::vector<char>
312  const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
313  const typename Triangulation<dim, spacedim>::CellStatus status)
314  {
315  typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
316  dof_handler);
317 
318  // create buffer for each individual object
319  std::vector<::Vector<typename VectorType::value_type>> dof_values(
320  input_vectors.size());
321 
322  unsigned int fe_index = 0;
323  if (dof_handler->has_hp_capabilities())
324  {
325  switch (status)
326  {
328  spacedim>::CELL_PERSIST:
330  spacedim>::CELL_REFINE:
331  {
332  fe_index = cell->future_fe_index();
333  break;
334  }
335 
337  spacedim>::CELL_COARSEN:
338  {
339  // In case of coarsening, we need to find a suitable FE index
340  // for the parent cell. We choose the 'least dominant fe'
341  // on all children from the associated FECollection.
342 # ifdef DEBUG
343  for (const auto &child : cell->child_iterators())
344  Assert(child->is_active() && child->coarsen_flag_set(),
345  typename ::Triangulation<
346  dim>::ExcInconsistentCoarseningFlags());
347 # endif
348 
349  fe_index = ::internal::hp::DoFHandlerImplementation::
350  dominated_future_fe_on_children<dim, spacedim>(cell);
351  break;
352  }
353 
354  default:
355  Assert(false, ExcInternalError());
356  break;
357  }
358  }
359 
360  const unsigned int dofs_per_cell =
361  dof_handler->get_fe(fe_index).n_dofs_per_cell();
362 
363  if (dofs_per_cell == 0)
364  return std::vector<char>(); // nothing to do for FE_Nothing
365 
366  auto it_input = input_vectors.cbegin();
367  auto it_output = dof_values.begin();
368  for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
369  {
370  it_output->reinit(dofs_per_cell);
371  cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
372  }
373 
374  return pack_dof_values<typename VectorType::value_type>(dof_values,
375  dofs_per_cell);
376  }
377 
378 
379 
380  template <int dim, typename VectorType, int spacedim>
381  void
383  const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
384  const typename Triangulation<dim, spacedim>::CellStatus status,
385  const boost::iterator_range<std::vector<char>::const_iterator>
386  & data_range,
387  std::vector<VectorType *> &all_out,
388  VectorType & valence)
389  {
390  typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
391  dof_handler);
392 
393  unsigned int fe_index = 0;
394  if (dof_handler->has_hp_capabilities())
395  {
396  switch (status)
397  {
399  spacedim>::CELL_PERSIST:
401  spacedim>::CELL_COARSEN:
402  {
403  fe_index = cell->active_fe_index();
404  break;
405  }
406 
408  spacedim>::CELL_REFINE:
409  {
410  // After refinement, this particular cell is no longer active,
411  // and its children have inherited its FE index. However, to
412  // unpack the data on the old cell, we need to recover its FE
413  // index from one of the children. Just to be sure, we also
414  // check if all children have the same FE index.
415  fe_index = cell->child(0)->active_fe_index();
416  for (unsigned int child_index = 1;
417  child_index < cell->n_children();
418  ++child_index)
419  Assert(cell->child(child_index)->active_fe_index() ==
420  fe_index,
421  ExcInternalError());
422  break;
423  }
424 
425  default:
426  Assert(false, ExcInternalError());
427  break;
428  }
429  }
430 
431  const unsigned int dofs_per_cell =
432  dof_handler->get_fe(fe_index).n_dofs_per_cell();
433 
434  if (dofs_per_cell == 0)
435  return; // nothing to do for FE_Nothing
436 
437  const std::vector<::Vector<typename VectorType::value_type>>
438  dof_values =
439  unpack_dof_values<typename VectorType::value_type>(data_range,
440  dofs_per_cell);
441 
442  // check if sizes match
443  Assert(dof_values.size() == all_out.size(), ExcInternalError());
444 
445  // check if we have enough dofs provided by the FE object
446  // to interpolate the transferred data correctly
447  for (auto it_dof_values = dof_values.begin();
448  it_dof_values != dof_values.end();
449  ++it_dof_values)
450  Assert(
451  dofs_per_cell == it_dof_values->size(),
452  ExcMessage(
453  "The transferred data was packed with a different number of dofs than the "
454  "currently registered FE object assigned to the DoFHandler has."));
455 
456  // distribute data for each registered vector on mesh
457  auto it_input = dof_values.cbegin();
458  auto it_output = all_out.begin();
459  for (; it_input != dof_values.cend(); ++it_input, ++it_output)
460  if (average_values)
461  cell->distribute_local_to_global_by_interpolation(*it_input,
462  *(*it_output),
463  fe_index);
464  else
465  cell->set_dof_values_by_interpolation(*it_input,
466  *(*it_output),
467  fe_index,
468  true);
469 
470  if (average_values)
471  {
472  // compute valence vector if averaging should be performed
473  Vector<typename VectorType::value_type> ones(dofs_per_cell);
474  ones = 1.0;
475  cell->distribute_local_to_global_by_interpolation(ones,
476  valence,
477  fe_index);
478  }
479  }
480  } // namespace distributed
481 } // namespace parallel
482 
483 
484 // explicit instantiations
485 # include "solution_transfer.inst"
486 
488 
489 #endif
virtual void clear()
Definition: vector.h:109
void unpack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range, std::vector< VectorType * > &all_out, VectorType &valence)
std::vector< char > pack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status)
SmartPointer< const DoFHandler< dim, spacedim >, SolutionTransfer< dim, VectorType, spacedim > > dof_handler
void prepare_for_coarsening_and_refinement(const std::vector< const VectorType * > &all_in)
void prepare_for_serialization(const VectorType &in)
void interpolate(std::vector< VectorType * > &all_out)
SolutionTransfer(const DoFHandler< dim, spacedim > &dof_handler, const bool average_values=false)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1509
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1703
static ::ExceptionBase & ExcMessage(std::string arg1)
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
unsigned short int fe_index
Definition: types.h:59
const ::Triangulation< dim, spacedim > & tria