Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10: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\}}\)
Loading...
Searching...
No Matches
solution_transfer.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2010 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#include <deal.II/base/config.h>
17
18#ifdef DEAL_II_WITH_P4EST
19
22
25
28
36# include <deal.II/lac/vector.h>
37
38# include <functional>
39# include <numeric>
40
41
43
44
45namespace
46{
55 template <typename value_type>
56 std::vector<char>
57 pack_dof_values(std::vector<Vector<value_type>> &dof_values,
58 const unsigned int dofs_per_cell)
59 {
60 for (const auto &values : dof_values)
61 {
62 AssertDimension(values.size(), dofs_per_cell);
63 (void)values;
64 }
65
66 const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
67
68 std::vector<char> buffer(dof_values.size() * bytes_per_entry);
69 for (unsigned int i = 0; i < dof_values.size(); ++i)
70 std::memcpy(&buffer[i * bytes_per_entry],
71 &dof_values[i](0),
72 bytes_per_entry);
73
74 return buffer;
75 }
76
77
78
82 template <typename value_type>
83 std::vector<Vector<value_type>>
84 unpack_dof_values(
85 const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
86 const unsigned int dofs_per_cell)
87 {
88 const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
89 const unsigned int n_elements = data_range.size() / bytes_per_entry;
90
91 Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
92
93 std::vector<Vector<value_type>> unpacked_data;
94 unpacked_data.reserve(n_elements);
95 for (unsigned int i = 0; i < n_elements; ++i)
96 {
97 Vector<value_type> dof_values(dofs_per_cell);
98 std::memcpy(&dof_values(0),
99 &(*std::next(data_range.begin(), i * bytes_per_entry)),
100 bytes_per_entry);
101 unpacked_data.emplace_back(std::move(dof_values));
102 }
103
104 return unpacked_data;
105 }
106} // namespace
107
108
109
110namespace parallel
111{
112 namespace distributed
113 {
114 template <int dim, typename VectorType, int spacedim>
116 const DoFHandler<dim, spacedim> &dof,
117 const bool average_values)
118 : dof_handler(&dof, typeid(*this).name())
119 , average_values(average_values)
120 , handle(numbers::invalid_unsigned_int)
121 {
122 Assert(
123 (dynamic_cast<
125 &dof_handler->get_triangulation()) != nullptr),
127 "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
128 }
129
130
131
132 template <int dim, typename VectorType, int spacedim>
133 void
136 const std::vector<const VectorType *> &all_in)
137 {
138 for (unsigned int i = 0; i < all_in.size(); ++i)
139 Assert(all_in[i]->size() == dof_handler->n_dofs(),
140 ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
141
142 input_vectors = all_in;
143 register_data_attach();
144 }
145
146
147
148 template <int dim, typename VectorType, int spacedim>
149 void
151 {
152 // TODO: casting away constness is bad
155 const_cast<::Triangulation<dim, spacedim> *>(
156 &dof_handler->get_triangulation())));
157 Assert(tria != nullptr, ExcInternalError());
158
160 ExcMessage("You can only add one solution per "
161 "SolutionTransfer object."));
162
163 handle = tria->register_data_attach(
164 [this](
166 const CellStatus status) {
167 return this->pack_callback(cell_, status);
168 },
169 /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
170 }
171
172
173
174 template <int dim, typename VectorType, int spacedim>
175 void
177 prepare_for_coarsening_and_refinement(const VectorType &in)
178 {
179 std::vector<const VectorType *> all_in(1, &in);
180 prepare_for_coarsening_and_refinement(all_in);
181 }
182
183
184
185 template <int dim, typename VectorType, int spacedim>
186 void
188 const VectorType &in)
189 {
190 std::vector<const VectorType *> all_in(1, &in);
191 prepare_for_serialization(all_in);
192 }
193
194
195
196 template <int dim, typename VectorType, int spacedim>
197 void
199 const std::vector<const VectorType *> &all_in)
200 {
201 prepare_for_coarsening_and_refinement(all_in);
202 }
203
204
205
206 template <int dim, typename VectorType, int spacedim>
207 void
209 {
210 std::vector<VectorType *> all_in(1, &in);
211 deserialize(all_in);
212 }
213
214
215
216 template <int dim, typename VectorType, int spacedim>
217 void
219 std::vector<VectorType *> &all_in)
220 {
221 register_data_attach();
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, int spacedim>
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
244 const_cast<::Triangulation<dim, spacedim> *>(
245 &dof_handler->get_triangulation())));
246 Assert(tria != nullptr, ExcInternalError());
247
248 if (average_values)
249 for (auto *const vec : all_out)
250 *vec = 0.0;
251
252 VectorType valence;
253
254 // initialize valence vector only if we need to average
255 if (average_values)
256 valence.reinit(*all_out[0]);
257
259 handle,
260 [this, &all_out, &valence](
262 const CellStatus status,
263 const boost::iterator_range<std::vector<char>::const_iterator>
264 &data_range) {
265 this->unpack_callback(cell_, status, data_range, all_out, valence);
266 });
267
268 if (average_values)
269 {
270 // finalize valence: compress and invert
271 using Number = typename VectorType::value_type;
272 valence.compress(VectorOperation::add);
273 for (const auto i : valence.locally_owned_elements())
274 valence[i] = (static_cast<Number>(valence[i]) == Number() ?
275 Number() :
276 (Number(1.0) / static_cast<Number>(valence[i])));
277 valence.compress(VectorOperation::insert);
278
279 for (auto *const vec : all_out)
280 {
281 // compress and weight with valence
282 vec->compress(VectorOperation::add);
283 vec->scale(valence);
284 }
285 }
286 else
287 {
288 for (auto *const vec : all_out)
289 vec->compress(VectorOperation::insert);
290 }
291
292 input_vectors.clear();
294 }
295
296
297
298 template <int dim, typename VectorType, int spacedim>
299 void
301 {
302 std::vector<VectorType *> all_out(1, &out);
303 interpolate(all_out);
304 }
305
306
307
308 template <int dim, typename VectorType, int spacedim>
309 std::vector<char>
312 const CellStatus status)
313 {
314 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
315 dof_handler);
316
317 // create buffer for each individual object
318 std::vector<::Vector<typename VectorType::value_type>> dof_values(
319 input_vectors.size());
320
321 unsigned int fe_index = 0;
322 if (dof_handler->has_hp_capabilities())
323 {
324 switch (status)
325 {
328 {
329 fe_index = cell->future_fe_index();
330 break;
331 }
332
334 {
335 // In case of coarsening, we need to find a suitable FE index
336 // for the parent cell. We choose the 'least dominant fe'
337 // on all children from the associated FECollection.
338# ifdef DEBUG
339 for (const auto &child : cell->child_iterators())
340 Assert(child->is_active() && child->coarsen_flag_set(),
341 typename ::Triangulation<
342 dim>::ExcInconsistentCoarseningFlags());
343# endif
344
345 fe_index = ::internal::hp::DoFHandlerImplementation::
346 dominated_future_fe_on_children<dim, spacedim>(cell);
347 break;
348 }
349
350 default:
352 break;
353 }
354 }
355
356 const unsigned int dofs_per_cell =
357 dof_handler->get_fe(fe_index).n_dofs_per_cell();
358
359 if (dofs_per_cell == 0)
360 return std::vector<char>(); // nothing to do for FE_Nothing
361
362 auto it_input = input_vectors.cbegin();
363 auto it_output = dof_values.begin();
364 for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
365 {
366 it_output->reinit(dofs_per_cell);
367 cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
368 }
369
370 return pack_dof_values<typename VectorType::value_type>(dof_values,
371 dofs_per_cell);
372 }
373
374
375
376 template <int dim, typename VectorType, int spacedim>
377 void
380 const CellStatus status,
381 const boost::iterator_range<std::vector<char>::const_iterator>
382 &data_range,
383 std::vector<VectorType *> &all_out,
384 VectorType &valence)
385 {
386 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
387 dof_handler);
388
389 unsigned int fe_index = 0;
390 if (dof_handler->has_hp_capabilities())
391 {
392 switch (status)
393 {
396 {
397 fe_index = cell->active_fe_index();
398 break;
399 }
400
402 {
403 // After refinement, this particular cell is no longer active,
404 // and its children have inherited its FE index. However, to
405 // unpack the data on the old cell, we need to recover its FE
406 // index from one of the children. Just to be sure, we also
407 // check if all children have the same FE index.
408 fe_index = cell->child(0)->active_fe_index();
409 for (unsigned int child_index = 1;
410 child_index < cell->n_children();
411 ++child_index)
412 Assert(cell->child(child_index)->active_fe_index() ==
413 fe_index,
415 break;
416 }
417
418 default:
420 break;
421 }
422 }
423
424 const unsigned int dofs_per_cell =
425 dof_handler->get_fe(fe_index).n_dofs_per_cell();
426
427 if (dofs_per_cell == 0)
428 return; // nothing to do for FE_Nothing
429
430 const std::vector<::Vector<typename VectorType::value_type>>
431 dof_values =
432 unpack_dof_values<typename VectorType::value_type>(data_range,
433 dofs_per_cell);
434
435 // check if sizes match
436 Assert(dof_values.size() == all_out.size(), ExcInternalError());
437
438 // check if we have enough dofs provided by the FE object
439 // to interpolate the transferred data correctly
440 for (auto it_dof_values = dof_values.begin();
441 it_dof_values != dof_values.end();
442 ++it_dof_values)
443 Assert(
444 dofs_per_cell == it_dof_values->size(),
446 "The transferred data was packed with a different number of dofs than the "
447 "currently registered FE object assigned to the DoFHandler has."));
448
449 // distribute data for each registered vector on mesh
450 auto it_input = dof_values.cbegin();
451 auto it_output = all_out.begin();
452 for (; it_input != dof_values.cend(); ++it_input, ++it_output)
453 if (average_values)
454 cell->distribute_local_to_global_by_interpolation(*it_input,
455 *(*it_output),
456 fe_index);
457 else
458 cell->set_dof_values_by_interpolation(*it_input,
459 *(*it_output),
460 fe_index,
461 true);
462
463 if (average_values)
464 {
465 // compute valence vector if averaging should be performed
466 Vector<typename VectorType::value_type> ones(dofs_per_cell);
467 ones = 1.0;
468 cell->distribute_local_to_global_by_interpolation(ones,
469 valence,
470 fe_index);
471 }
472 }
473 } // namespace distributed
474} // namespace parallel
475
476
477// explicit instantiations
478# include "solution_transfer.inst"
479
481
482#endif
CellStatus
Definition cell_status.h:31
@ cell_will_be_refined
@ children_will_be_coarsened
unsigned int register_data_attach(const std::function< std::vector< char >(const cell_iterator &, const ::CellStatus)> &pack_callback, const bool returns_variable_size_data)
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)
std::vector< char > pack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellStatus status)
void unpack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range, std::vector< VectorType * > &all_out, VectorType &valence)
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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:269
#define DEAL_II_ASSERT_UNREACHABLE()
static const unsigned int invalid_unsigned_int
Definition types.h:220
const InputIterator OutputIterator out
Definition parallel.h:167