Reference documentation for deal.II version 9.4.0
\(\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
26
29
31
39# include <deal.II/lac/vector.h>
40
41# include <functional>
42# include <numeric>
43
44
46
47
48namespace
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
113namespace parallel
114{
115 namespace distributed
116 {
117 template <int dim, typename VectorType, int spacedim>
119 const DoFHandler<dim, spacedim> &dof,
120 const bool average_values)
121 : dof_handler(&dof, typeid(*this).name())
122 , average_values(average_values)
124 {
125 Assert(
126 (dynamic_cast<
128 &dof_handler->get_triangulation()) != nullptr),
130 "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
131 }
132
133
134
135 template <int dim, typename VectorType, int spacedim>
136 void
139 const std::vector<const VectorType *> &all_in)
140 {
141 for (unsigned int i = 0; i < all_in.size(); ++i)
142 Assert(all_in[i]->size() == dof_handler->n_dofs(),
143 ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
144
145 input_vectors = all_in;
146 register_data_attach();
147 }
148
149
150
151 template <int dim, typename VectorType, int spacedim>
152 void
154 {
155 // TODO: casting away constness is bad
158 const_cast<::Triangulation<dim, spacedim> *>(
159 &dof_handler->get_triangulation())));
160 Assert(tria != nullptr, ExcInternalError());
161
162 handle = tria->register_data_attach(
163 [this](
165 const typename Triangulation<dim, spacedim>::CellStatus status) {
166 return this->pack_callback(cell_, status);
167 },
168 /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
169 }
170
171
172
173 template <int dim, typename VectorType, int spacedim>
174 void
176 prepare_for_coarsening_and_refinement(const VectorType &in)
177 {
178 std::vector<const VectorType *> all_in(1, &in);
179 prepare_for_coarsening_and_refinement(all_in);
180 }
181
182
183
184 template <int dim, typename VectorType, int spacedim>
185 void
187 const VectorType &in)
188 {
189 std::vector<const VectorType *> all_in(1, &in);
190 prepare_for_serialization(all_in);
191 }
192
193
194
195 template <int dim, typename VectorType, int spacedim>
196 void
198 const std::vector<const VectorType *> &all_in)
199 {
200 prepare_for_coarsening_and_refinement(all_in);
201 }
202
203
204
205 template <int dim, typename VectorType, int spacedim>
206 void
208 {
209 std::vector<VectorType *> all_in(1, &in);
210 deserialize(all_in);
211 }
212
213
214
215 template <int dim, typename VectorType, int spacedim>
216 void
218 std::vector<VectorType *> &all_in)
219 {
220 register_data_attach();
221
222 // this makes interpolate() happy
223 input_vectors.resize(all_in.size());
224
225 interpolate(all_in);
226 }
227
228
229 template <int dim, typename VectorType, int spacedim>
230 void
232 std::vector<VectorType *> &all_out)
233 {
234 Assert(input_vectors.size() == all_out.size(),
235 ExcDimensionMismatch(input_vectors.size(), all_out.size()));
236 for (unsigned int i = 0; i < all_out.size(); ++i)
237 Assert(all_out[i]->size() == dof_handler->n_dofs(),
238 ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
239
240 // TODO: casting away constness is bad
243 const_cast<::Triangulation<dim, spacedim> *>(
244 &dof_handler->get_triangulation())));
245 Assert(tria != nullptr, ExcInternalError());
246
247 if (average_values)
248 for (const auto vec : all_out)
249 *vec = 0.0;
250
251 VectorType valence;
252
253 // initialize valence vector only if we need to average
254 if (average_values)
255 valence.reinit(*all_out[0]);
256
257 tria->notify_ready_to_unpack(
258 handle,
259 [this, &all_out, &valence](
261 const typename Triangulation<dim, spacedim>::CellStatus status,
262 const boost::iterator_range<std::vector<char>::const_iterator>
263 &data_range) {
264 this->unpack_callback(cell_, status, data_range, all_out, valence);
265 });
266
267 if (average_values)
268 {
269 // finalize valence: compress and invert
270 using Number = typename VectorType::value_type;
271 valence.compress(::VectorOperation::add);
272 for (const auto i : valence.locally_owned_elements())
273 valence[i] = (static_cast<Number>(valence[i]) == Number() ?
274 Number() :
275 (Number(1.0) / static_cast<Number>(valence[i])));
276 valence.compress(::VectorOperation::insert);
277
278 for (const auto vec : all_out)
279 {
280 // compress and weight with valence
281 vec->compress(::VectorOperation::add);
282 vec->scale(valence);
283 }
284 }
285 else
286 {
287 for (const auto vec : all_out)
288 vec->compress(::VectorOperation::insert);
289 }
290
291 input_vectors.clear();
292 }
293
294
295
296 template <int dim, typename VectorType, int spacedim>
297 void
299 {
300 std::vector<VectorType *> all_out(1, &out);
301 interpolate(all_out);
302 }
303
304
305
306 template <int dim, typename VectorType, int spacedim>
307 std::vector<char>
310 const typename Triangulation<dim, spacedim>::CellStatus status)
311 {
312 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
313 dof_handler);
314
315 // create buffer for each individual object
316 std::vector<::Vector<typename VectorType::value_type>> dof_values(
317 input_vectors.size());
318
319 unsigned int fe_index = 0;
320 if (dof_handler->has_hp_capabilities())
321 {
322 switch (status)
323 {
325 spacedim>::CELL_PERSIST:
327 spacedim>::CELL_REFINE:
328 {
329 fe_index = cell->future_fe_index();
330 break;
331 }
332
334 spacedim>::CELL_COARSEN:
335 {
336 // In case of coarsening, we need to find a suitable FE index
337 // for the parent cell. We choose the 'least dominant fe'
338 // on all children from the associated FECollection.
339# ifdef DEBUG
340 for (const auto &child : cell->child_iterators())
341 Assert(child->is_active() && child->coarsen_flag_set(),
342 typename ::Triangulation<
343 dim>::ExcInconsistentCoarseningFlags());
344# endif
345
346 fe_index = ::internal::hp::DoFHandlerImplementation::
347 dominated_future_fe_on_children<dim, spacedim>(cell);
348 break;
349 }
350
351 default:
352 Assert(false, ExcInternalError());
353 break;
354 }
355 }
356
357 const unsigned int dofs_per_cell =
358 dof_handler->get_fe(fe_index).n_dofs_per_cell();
359
360 if (dofs_per_cell == 0)
361 return std::vector<char>(); // nothing to do for FE_Nothing
362
363 auto it_input = input_vectors.cbegin();
364 auto it_output = dof_values.begin();
365 for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
366 {
367 it_output->reinit(dofs_per_cell);
368 cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
369 }
370
371 return pack_dof_values<typename VectorType::value_type>(dof_values,
372 dofs_per_cell);
373 }
374
375
376
377 template <int dim, typename VectorType, int spacedim>
378 void
381 const typename Triangulation<dim, spacedim>::CellStatus status,
382 const boost::iterator_range<std::vector<char>::const_iterator>
383 & data_range,
384 std::vector<VectorType *> &all_out,
385 VectorType & valence)
386 {
387 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
388 dof_handler);
389
390 unsigned int fe_index = 0;
391 if (dof_handler->has_hp_capabilities())
392 {
393 switch (status)
394 {
396 spacedim>::CELL_PERSIST:
398 spacedim>::CELL_COARSEN:
399 {
400 fe_index = cell->active_fe_index();
401 break;
402 }
403
405 spacedim>::CELL_REFINE:
406 {
407 // After refinement, this particular cell is no longer active,
408 // and its children have inherited its FE index. However, to
409 // unpack the data on the old cell, we need to recover its FE
410 // index from one of the children. Just to be sure, we also
411 // check if all children have the same FE index.
412 fe_index = cell->child(0)->active_fe_index();
413 for (unsigned int child_index = 1;
414 child_index < cell->n_children();
415 ++child_index)
416 Assert(cell->child(child_index)->active_fe_index() ==
417 fe_index,
419 break;
420 }
421
422 default:
423 Assert(false, ExcInternalError());
424 break;
425 }
426 }
427
428 const unsigned int dofs_per_cell =
429 dof_handler->get_fe(fe_index).n_dofs_per_cell();
430
431 if (dofs_per_cell == 0)
432 return; // nothing to do for FE_Nothing
433
434 const std::vector<::Vector<typename VectorType::value_type>>
435 dof_values =
436 unpack_dof_values<typename VectorType::value_type>(data_range,
437 dofs_per_cell);
438
439 // check if sizes match
440 Assert(dof_values.size() == all_out.size(), ExcInternalError());
441
442 // check if we have enough dofs provided by the FE object
443 // to interpolate the transferred data correctly
444 for (auto it_dof_values = dof_values.begin();
445 it_dof_values != dof_values.end();
446 ++it_dof_values)
447 Assert(
448 dofs_per_cell == it_dof_values->size(),
450 "The transferred data was packed with a different number of dofs than the "
451 "currently registered FE object assigned to the DoFHandler has."));
452
453 // distribute data for each registered vector on mesh
454 auto it_input = dof_values.cbegin();
455 auto it_output = all_out.begin();
456 for (; it_input != dof_values.cend(); ++it_input, ++it_output)
457 if (average_values)
458 cell->distribute_local_to_global_by_interpolation(*it_input,
459 *(*it_output),
460 fe_index);
461 else
462 cell->set_dof_values_by_interpolation(*it_input,
463 *(*it_output),
464 fe_index,
465 true);
466
467 if (average_values)
468 {
469 // compute valence vector if averaging should be performed
470 Vector<typename VectorType::value_type> ones(dofs_per_cell);
471 ones = 1.0;
472 cell->distribute_local_to_global_by_interpolation(ones,
473 valence,
474 fe_index);
475 }
476 }
477 } // namespace distributed
478} // namespace parallel
479
480
481// explicit instantiations
482# include "solution_transfer.inst"
483
485
486#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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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:201
const ::Triangulation< dim, spacedim > & tria