Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
solution_transfer.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2009 - 2023 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
37# include <deal.II/lac/vector.h>
38
39# include <functional>
40# include <numeric>
41
42
44
45
46namespace
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
111namespace 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),
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](
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](
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>
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
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,
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(),
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()
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)
typename ::Triangulation< dim, spacedim >::CellStatus CellStatus
Definition tria.h:294
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition tria.h:270
static const unsigned int invalid_unsigned_int
Definition types.h:213
const ::Triangulation< dim, spacedim > & tria