Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2968-g5f01c80b02 2025-03-29 13:10:00+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
solution_transfer.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 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#include <deal.II/base/config.h>
16
18
21
25
26#include <deal.II/fe/fe.h>
27
28#include <deal.II/grid/tria.h>
31
39#include <deal.II/lac/vector.h>
41
43
44#include <boost/range/iterator_range_core.hpp>
45
46#include <functional>
47#include <numeric>
48
49
51
52
53namespace
54{
63 template <typename value_type>
64 std::vector<char>
65 pack_dof_values(std::vector<Vector<value_type>> &dof_values,
66 const unsigned int dofs_per_cell)
67 {
68 for (const auto &values : dof_values)
69 {
70 AssertDimension(values.size(), dofs_per_cell);
71 (void)values;
72 }
73
74 const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
75
76 std::vector<char> buffer(dof_values.size() * bytes_per_entry);
77 for (unsigned int i = 0; i < dof_values.size(); ++i)
78 std::memcpy(&buffer[i * bytes_per_entry],
79 &dof_values[i](0),
80 bytes_per_entry);
81
82 return buffer;
83 }
84
85
86
90 template <typename value_type>
91 std::vector<Vector<value_type>>
92 unpack_dof_values(
93 const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
94 const unsigned int dofs_per_cell)
95 {
96 const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
97 const unsigned int n_elements = data_range.size() / bytes_per_entry;
98
99 Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
100
101 std::vector<Vector<value_type>> unpacked_data;
102 unpacked_data.reserve(n_elements);
103 for (unsigned int i = 0; i < n_elements; ++i)
104 {
105 Vector<value_type> dof_values(dofs_per_cell);
106 std::memcpy(&dof_values(0),
107 &(*std::next(data_range.begin(), i * bytes_per_entry)),
108 bytes_per_entry);
109 unpacked_data.emplace_back(std::move(dof_values));
110 }
111
112 return unpacked_data;
113 }
114} // namespace
115
116
117
118template <int dim, typename VectorType, int spacedim>
120 const DoFHandler<dim, spacedim> &dof,
121 const bool average_values)
122 : dof_handler(&dof, typeid(*this).name())
123 , average_values(average_values)
124 , handle(numbers::invalid_unsigned_int)
125{}
126
127
128
129template <int dim, typename VectorType, int spacedim>
130void
133 const std::vector<const VectorType *> &all_in)
134{
135 const ::internal::parallel::shared::
136 TemporarilyRestoreSubdomainIds<dim, spacedim>
137 subdomain_modifier(dof_handler->get_triangulation());
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
149template <int dim, typename VectorType, int spacedim>
150void
152 prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
153{
154 std::vector<const VectorType *> temp(all_in.size());
155
156 for (std::size_t i = 0; i < temp.size(); ++i)
157 temp[i] = &(all_in[i]);
158
159 this->prepare_for_coarsening_and_refinement(temp);
160}
161
162
163
164template <int dim, typename VectorType, int spacedim>
165void
167{
168 // TODO: casting away constness is bad
169 auto tria = const_cast<::Triangulation<dim, spacedim> *>(
170 &dof_handler->get_triangulation());
171 Assert(tria != nullptr, ExcInternalError());
172
174 ExcMessage("You can only add one solution per "
175 "SolutionTransfer object."));
176
177 handle = tria->register_data_attach(
178 [this](const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
179 const CellStatus status) {
180 return this->pack_callback(cell_, status);
181 },
182 /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
183}
184
185
186
187template <int dim, typename VectorType, int spacedim>
188void
190 prepare_for_coarsening_and_refinement(const VectorType &in)
191{
192 std::vector<const VectorType *> all_in(1, &in);
193 prepare_for_coarsening_and_refinement(all_in);
194}
195
196
197
198template <int dim, typename VectorType, int spacedim>
199void
201 const VectorType &in)
202{
203 std::vector<const VectorType *> all_in(1, &in);
204 prepare_for_serialization(all_in);
205}
206
207
208
209template <int dim, typename VectorType, int spacedim>
210void
212 const std::vector<const VectorType *> &all_in)
213{
214 prepare_for_coarsening_and_refinement(all_in);
215}
216
217
218
219template <int dim, typename VectorType, int spacedim>
220void
222{
223 std::vector<VectorType *> all_in(1, &in);
224 deserialize(all_in);
225}
226
227
228
229template <int dim, typename VectorType, int spacedim>
230void
232 std::vector<VectorType *> &all_in)
233{
234 register_data_attach();
235
236 // this makes interpolate() happy
237 input_vectors.resize(all_in.size());
238
239 interpolate(all_in);
240}
241
242
243template <int dim, typename VectorType, int spacedim>
244void
246 std::vector<VectorType *> &all_out)
247{
248 const ::internal::parallel::shared::
249 TemporarilyRestoreSubdomainIds<dim, spacedim>
250 subdomain_modifier(dof_handler->get_triangulation());
251
252 Assert(input_vectors.size() == all_out.size(),
253 ExcDimensionMismatch(input_vectors.size(), all_out.size()));
254 for (unsigned int i = 0; i < all_out.size(); ++i)
255 Assert(all_out[i]->size() == dof_handler->n_dofs(),
256 ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
257
258 // TODO: casting away constness is bad
259 auto tria = const_cast<::Triangulation<dim, spacedim> *>(
260 &dof_handler->get_triangulation());
261 Assert(tria != nullptr, ExcInternalError());
262 Assert(
265 "You can only call interpolate() once per SolutionTransfer object."));
266
267 using Number = typename VectorType::value_type;
268
269 if (average_values)
270 for (auto *const vec : all_out)
271 *vec = Number();
272
273 VectorType valence;
274
275 // initialize valence vector only if we need to average
276 if (average_values)
277 valence.reinit(*all_out[0]);
278
279 tria->notify_ready_to_unpack(
280 handle,
281 [this, &all_out, &valence](
283 const CellStatus status,
284 const boost::iterator_range<std::vector<char>::const_iterator>
285 &data_range) {
286 this->unpack_callback(cell_, status, data_range, all_out, valence);
287 });
288
289 if (average_values)
290 {
291 // finalize valence: compress and invert
292 valence.compress(VectorOperation::add);
293 for (const auto i : valence.locally_owned_elements())
294 valence[i] = (static_cast<Number>(valence[i]) == Number() ?
295 Number() :
296 (Number(1.0) / static_cast<Number>(valence[i])));
297 valence.compress(VectorOperation::insert);
298
299 for (auto *const vec : all_out)
300 {
301 // compress and weight with valence
302 vec->compress(VectorOperation::add);
303 vec->scale(valence);
304 }
305 }
306 else
307 {
308 for (auto *const vec : all_out)
309 vec->compress(VectorOperation::insert);
310 }
311
312 input_vectors.clear();
314}
315
316
317template <int dim, typename VectorType, int spacedim>
318void
320 std::vector<VectorType> &all_out)
321{
322 std::vector<VectorType *> temp(all_out.size());
323
324 for (std::size_t i = 0; i < temp.size(); ++i)
325 temp[i] = &(all_out[i]);
326
327 this->interpolate(temp);
328}
329
330
331
332template <int dim, typename VectorType, int spacedim>
333void
335{
336 std::vector<VectorType *> all_out(1, &out);
337 interpolate(all_out);
338}
339
340
341
342template <int dim, typename VectorType, int spacedim>
343std::vector<char>
346 const CellStatus status)
347{
348 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
349
350 // create buffer for each individual object
351 std::vector<::Vector<typename VectorType::value_type>> dof_values(
352 input_vectors.size());
353
354 unsigned int fe_index = 0;
355 if (dof_handler->has_hp_capabilities())
356 {
357 switch (status)
358 {
361 {
362 fe_index = cell->future_fe_index();
363 break;
364 }
365
367 {
368 // In case of coarsening, we need to find a suitable FE index
369 // for the parent cell. We choose the 'least dominant fe'
370 // on all children from the associated FECollection.
371 if constexpr (running_in_debug_mode())
372 {
373 for (const auto &child : cell->child_iterators())
374 Assert(child->is_active() && child->coarsen_flag_set(),
375 typename ::Triangulation<
376 dim>::ExcInconsistentCoarseningFlags());
377 }
378
379 fe_index = ::internal::hp::DoFHandlerImplementation::
380 dominated_future_fe_on_children<dim, spacedim>(cell);
381 break;
382 }
383
384 default:
385 Assert(false, ExcInternalError());
386 break;
387 }
388 }
389
390 const unsigned int dofs_per_cell =
391 dof_handler->get_fe(fe_index).n_dofs_per_cell();
392
393 if (dofs_per_cell == 0)
394 return std::vector<char>(); // nothing to do for FE_Nothing
395
396 auto it_input = input_vectors.cbegin();
397 auto it_output = dof_values.begin();
398 for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
399 {
400 it_output->reinit(dofs_per_cell);
401 cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
402 }
403
404 return pack_dof_values<typename VectorType::value_type>(dof_values,
405 dofs_per_cell);
406}
407
408
409
410template <int dim, typename VectorType, int spacedim>
411void
414 const CellStatus status,
415 const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
416 std::vector<VectorType *> &all_out,
417 VectorType &valence)
418{
419 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
420
421 unsigned int fe_index = 0;
422 if (dof_handler->has_hp_capabilities())
423 {
424 switch (status)
425 {
428 {
429 fe_index = cell->active_fe_index();
430 break;
431 }
432
434 {
435 // After refinement, this particular cell is no longer active,
436 // and its children have inherited its FE index. However, to
437 // unpack the data on the old cell, we need to recover its FE
438 // index from one of the children. Just to be sure, we also
439 // check if all children have the same FE index.
440 fe_index = cell->child(0)->active_fe_index();
441 for (unsigned int child_index = 1;
442 child_index < cell->n_children();
443 ++child_index)
444 Assert(cell->child(child_index)->active_fe_index() == fe_index,
446 break;
447 }
448
449 default:
450 Assert(false, ExcInternalError());
451 break;
452 }
453 }
454
455 const unsigned int dofs_per_cell =
456 dof_handler->get_fe(fe_index).n_dofs_per_cell();
457
458 if (dofs_per_cell == 0)
459 return; // nothing to do for FE_Nothing
460
461 const std::vector<::Vector<typename VectorType::value_type>>
462 dof_values =
463 unpack_dof_values<typename VectorType::value_type>(data_range,
464 dofs_per_cell);
465
466 // check if sizes match
467 AssertDimension(dof_values.size(), all_out.size());
468
469 // check if we have enough dofs provided by the FE object
470 // to interpolate the transferred data correctly
471 for (auto it_dof_values = dof_values.begin();
472 it_dof_values != dof_values.end();
473 ++it_dof_values)
474 Assert(
475 dofs_per_cell == it_dof_values->size(),
477 "The transferred data was packed with a different number of dofs than the "
478 "currently registered FE object assigned to the DoFHandler has."));
479
480 // distribute data for each registered vector on mesh
481 auto it_input = dof_values.cbegin();
482 auto it_output = all_out.begin();
483 for (; it_input != dof_values.cend(); ++it_input, ++it_output)
484 if (average_values)
485 cell->distribute_local_to_global_by_interpolation(*it_input,
486 *(*it_output),
487 fe_index);
488 else
489 cell->set_dof_values_by_interpolation(*it_input,
490 *(*it_output),
491 fe_index,
492 true);
493
494 if (average_values)
495 {
496 // compute valence vector if averaging should be performed
497 Vector<typename VectorType::value_type> ones(dofs_per_cell);
498 ones = 1.0;
499 cell->distribute_local_to_global_by_interpolation(ones,
500 valence,
501 fe_index);
502 }
503}
504
505
506
507template <int dim, typename VectorType, int spacedim>
508void
513
514
515
516namespace Legacy
517{
518
519 template <int dim, typename VectorType, int spacedim>
521 const DoFHandler<dim, spacedim> &dof)
522 : dof_handler(&dof, typeid(*this).name())
523 , n_dofs_old(0)
524 , prepared_for(none)
525 {
526 Assert((dynamic_cast<
528 &dof_handler->get_triangulation()) == nullptr),
530 "You are calling the ::SolutionTransfer class "
531 "with a DoFHandler that is built on a "
532 "parallel::distributed::Triangulation. This will not "
533 "work for parallel computations. You probably want to "
534 "use the parallel::distributed::SolutionTransfer class."));
535 }
536
537
538
539 template <int dim, typename VectorType, int spacedim>
544
545
546
547 template <int dim, typename VectorType, int spacedim>
548 void
550 {
551 indices_on_cell.clear();
552 dof_values_on_cell.clear();
553 cell_map.clear();
554
555 prepared_for = none;
556 }
557
558
559
560 template <int dim, typename VectorType, int spacedim>
561 void
563 {
564 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
565 Assert(prepared_for != coarsening_and_refinement,
566 ExcAlreadyPrepForCoarseAndRef());
567
568 clear();
569
570 // We need to access dof indices on the entire domain. For
571 // parallel::shared::Triangulations, ownership of cells might change. If
572 // they allow artificial cells, we need to restore the "true" cell owners
573 // temporarily.
574 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
575 // the current set of subdomain ids, set subdomain ids to the "true" owner
576 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
577 // object, and later restore these flags when it is destroyed.
578 const ::internal::parallel::shared::
579 TemporarilyRestoreSubdomainIds<dim, spacedim>
580 subdomain_modifier(dof_handler->get_triangulation());
581
582 const unsigned int n_active_cells =
583 dof_handler->get_triangulation().n_active_cells();
584 n_dofs_old = dof_handler->n_dofs();
585
586 // efficient reallocation of indices_on_cell
587 std::vector<std::vector<types::global_dof_index>>(n_active_cells)
588 .swap(indices_on_cell);
589
590 for (const auto &cell : dof_handler->active_cell_iterators())
591 {
592 const unsigned int i = cell->active_cell_index();
593 indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
594 // on each cell store the indices of the
595 // dofs. after refining we get the values
596 // on the children by taking these
597 // indices, getting the respective values
598 // out of the data vectors and prolonging
599 // them to the children
600 cell->get_dof_indices(indices_on_cell[i]);
601 cell_map[std::make_pair(cell->level(), cell->index())] =
602 Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
603 }
604 prepared_for = pure_refinement;
605 }
606
607
608
609 template <int dim, typename VectorType, int spacedim>
610 void
612 const VectorType &in,
613 VectorType &out) const
614 {
615 Assert(prepared_for == pure_refinement, ExcNotPrepared());
616 Assert(in.size() == n_dofs_old,
617 ExcDimensionMismatch(in.size(), n_dofs_old));
618 Assert(out.size() == dof_handler->n_dofs(),
619 ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
620 Assert(&in != &out,
621 ExcMessage("Vectors cannot be used as input and output"
622 " at the same time!"));
623
624 // We need to access dof indices on the entire domain. For
625 // parallel::shared::Triangulations, ownership of cells might change. If
626 // they allow artificial cells, we need to restore the "true" cell owners
627 // temporarily.
628 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
629 // the current set of subdomain ids, set subdomain ids to the "true" owner
630 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
631 // object, and later restore these flags when it is destroyed.
632 const ::internal::parallel::shared::
633 TemporarilyRestoreSubdomainIds<dim, spacedim>
634 subdomain_modifier(dof_handler->get_triangulation());
635
637
638 typename std::map<std::pair<unsigned int, unsigned int>,
639 Pointerstruct>::const_iterator pointerstruct,
640 cell_map_end = cell_map.end();
641
642 for (const auto &cell : dof_handler->cell_iterators())
643 {
644 pointerstruct =
645 cell_map.find(std::make_pair(cell->level(), cell->index()));
646
647 if (pointerstruct != cell_map_end)
648 // this cell was refined or not
649 // touched at all, so we can get
650 // the new values by just setting
651 // or interpolating to the children,
652 // which is both done by one
653 // function
654 {
655 const unsigned int this_fe_index =
656 pointerstruct->second.active_fe_index;
657 const unsigned int dofs_per_cell =
658 cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
659 local_values.reinit(dofs_per_cell, true);
660
661 // make sure that the size of the stored indices is the same as
662 // dofs_per_cell. since we store the desired fe_index, we know
663 // what this size should be
664 Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
666 for (unsigned int i = 0; i < dofs_per_cell; ++i)
667 local_values(i) =
669 in, (*pointerstruct->second.indices_ptr)[i]);
670 cell->set_dof_values_by_interpolation(local_values,
671 out,
672 this_fe_index,
673 true);
674 }
675 }
676 }
677
678
679
680 namespace internal
681 {
695 template <int dim, int spacedim>
696 void
698 const DoFHandler<dim, spacedim> &dof,
699 ::Table<2, FullMatrix<double>> &matrices)
700 {
701 if (dof.has_hp_capabilities() == false)
702 return;
703
704 const ::hp::FECollection<dim, spacedim> &fe =
705 dof.get_fe_collection();
706 matrices.reinit(fe.size(), fe.size());
707 for (unsigned int i = 0; i < fe.size(); ++i)
708 for (unsigned int j = 0; j < fe.size(); ++j)
709 if (i != j)
710 {
711 matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
712 fe[j].n_dofs_per_cell());
713
714 // see if we can get the interpolation matrices for this
715 // combination of elements. if not, reset the matrix sizes to zero
716 // to indicate that this particular combination isn't
717 // supported. this isn't an outright error right away since we may
718 // never need to actually interpolate between these two elements
719 // on actual cells; we simply have to trigger an error if someone
720 // actually tries
721 try
722 {
723 fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
724 }
725 catch (const typename FiniteElement<dim, spacedim>::
726 ExcInterpolationNotImplemented &)
727 {
728 matrices(i, j).reinit(0, 0);
729 }
730 }
731 }
732
733
734 template <int dim, int spacedim>
735 void
737 std::vector<std::vector<bool>> &)
738 {}
739
740 template <int dim, int spacedim>
741 void
743 const ::hp::FECollection<dim, spacedim> &fe,
744 std::vector<std::vector<bool>> &restriction_is_additive)
745 {
746 restriction_is_additive.resize(fe.size());
747 for (unsigned int f = 0; f < fe.size(); ++f)
748 {
749 restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
750 for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
751 restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
752 }
753 }
754 } // namespace internal
755
756
757
758 template <int dim, typename VectorType, int spacedim>
759 void
761 prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
762 {
763 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
764 Assert(prepared_for != coarsening_and_refinement,
765 ExcAlreadyPrepForCoarseAndRef());
766
767 clear();
768 n_dofs_old = dof_handler->n_dofs();
769 const unsigned int in_size = all_in.size();
770
771 if constexpr (running_in_debug_mode())
772 {
773 Assert(in_size != 0,
774 ExcMessage("The array of input vectors you pass to this "
775 "function has no elements. This is not useful."));
776 for (unsigned int i = 0; i < in_size; ++i)
777 {
778 Assert(all_in[i].size() == n_dofs_old,
779 ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
780 }
781 }
782
783 // We need to access dof indices on the entire domain. For
784 // parallel::shared::Triangulations, ownership of cells might change. If
785 // they allow artificial cells, we need to restore the "true" cell owners
786 // temporarily.
787 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
788 // the current set of subdomain ids, set subdomain ids to the "true" owner
789 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
790 // object, and later restore these flags when it is destroyed.
791 const ::internal::parallel::shared::
792 TemporarilyRestoreSubdomainIds<dim, spacedim>
793 subdomain_modifier(dof_handler->get_triangulation());
794
795 // first count the number
796 // of cells that will be coarsened
797 // and that'll stay or be refined
798 unsigned int n_cells_to_coarsen = 0;
799 unsigned int n_cells_to_stay_or_refine = 0;
800 for (const auto &act_cell : dof_handler->active_cell_iterators())
801 {
802 if (act_cell->coarsen_flag_set())
803 ++n_cells_to_coarsen;
804 else
805 ++n_cells_to_stay_or_refine;
806 }
807 Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
808 dof_handler->get_triangulation().n_active_cells(),
810
811 unsigned int n_coarsen_fathers = 0;
812 for (const auto &cell : dof_handler->cell_iterators())
813 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
814 ++n_coarsen_fathers;
815 Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
816 (void)n_cells_to_coarsen;
817
818 // allocate the needed memory. initialize
819 // the following arrays in an efficient
820 // way, without copying much
821 std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
822 .swap(indices_on_cell);
823
824 std::vector<std::vector<Vector<typename VectorType::value_type>>>(
825 n_coarsen_fathers,
826 std::vector<Vector<typename VectorType::value_type>>(in_size))
827 .swap(dof_values_on_cell);
828
829 Table<2, FullMatrix<double>> interpolation_hp;
830 std::vector<std::vector<bool>> restriction_is_additive;
831
832 internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
833 internal::restriction_additive(dof_handler->get_fe_collection(),
834 restriction_is_additive);
835
836 // we need counters for
837 // the 'to_stay_or_refine' cells 'n_sr' and
838 // the 'coarsen_fathers' cells 'n_cf',
839 unsigned int n_sr = 0, n_cf = 0;
840 for (const auto &cell : dof_handler->cell_iterators())
841 {
842 // CASE 1: active cell that remains as it is
843 if (cell->is_active() && !cell->coarsen_flag_set())
844 {
845 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
846 indices_on_cell[n_sr].resize(dofs_per_cell);
847 // cell will not be coarsened,
848 // so we get away by storing the
849 // dof indices and later
850 // interpolating to the children
851 cell->get_dof_indices(indices_on_cell[n_sr]);
852 cell_map[std::make_pair(cell->level(), cell->index())] =
853 Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
854 ++n_sr;
855 }
856
857 // CASE 2: cell is inactive but will become active
858 else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
859 {
860 // we will need to interpolate from the children of this cell
861 // to the current one. in the hp-context, this also means
862 // we need to figure out which finite element space to interpolate
863 // to since that is not implied by the global FE as in the non-hp-
864 // case. we choose the 'least dominant fe' on all children from
865 // the associated FECollection.
866 std::set<unsigned int> fe_indices_children;
867 for (const auto &child : cell->child_iterators())
868 {
869 Assert(child->is_active() && child->coarsen_flag_set(),
870 typename ::Triangulation<
871 dim>::ExcInconsistentCoarseningFlags());
872
873 fe_indices_children.insert(child->active_fe_index());
874 }
875 Assert(!fe_indices_children.empty(), ExcInternalError());
876
877 const unsigned int target_fe_index =
878 dof_handler->get_fe_collection().find_dominated_fe_extended(
879 fe_indices_children, /*codim=*/0);
880
881 Assert(target_fe_index != numbers::invalid_unsigned_int,
882 ::internal::hp::DoFHandlerImplementation::
883 ExcNoDominatedFiniteElementOnChildren());
884
885 const unsigned int dofs_per_cell =
886 dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
887
888 std::vector<Vector<typename VectorType::value_type>>(
889 in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
890 .swap(dof_values_on_cell[n_cf]);
891
892
893 // store the data of each of the input vectors. get this data
894 // as interpolated onto a finite element space that encompasses
895 // that of all the children. note that
896 // cell->get_interpolated_dof_values already does all of the
897 // interpolations between spaces
898 for (unsigned int j = 0; j < in_size; ++j)
899 cell->get_interpolated_dof_values(all_in[j],
900 dof_values_on_cell[n_cf][j],
901 target_fe_index);
902 cell_map[std::make_pair(cell->level(), cell->index())] =
903 Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
904 ++n_cf;
905 }
906 }
907 Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
908 Assert(n_cf == n_coarsen_fathers, ExcInternalError());
909
910 prepared_for = coarsening_and_refinement;
911 }
912
913
914
915 template <int dim, typename VectorType, int spacedim>
916 void
918 prepare_for_coarsening_and_refinement(const VectorType &in)
919 {
920 std::vector<VectorType> all_in(1, in);
921 prepare_for_coarsening_and_refinement(all_in);
922 }
923
924
925
926 template <int dim, typename VectorType, int spacedim>
927 void
929 const std::vector<VectorType> &all_in,
930 std::vector<VectorType> &all_out) const
931 {
932 const unsigned int size = all_in.size();
933 if constexpr (running_in_debug_mode())
934 {
935 Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
936 Assert(all_out.size() == size,
937 ExcDimensionMismatch(all_out.size(), size));
938 for (unsigned int i = 0; i < size; ++i)
939 Assert(all_in[i].size() == n_dofs_old,
940 ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
941 for (unsigned int i = 0; i < all_out.size(); ++i)
942 Assert(all_out[i].size() == dof_handler->n_dofs(),
943 ExcDimensionMismatch(all_out[i].size(),
944 dof_handler->n_dofs()));
945 for (unsigned int i = 0; i < size; ++i)
946 for (unsigned int j = 0; j < size; ++j)
947 Assert(&all_in[i] != &all_out[j],
948 ExcMessage("Vectors cannot be used as input and output"
949 " at the same time!"));
950 }
951
952 // We need to access dof indices on the entire domain. For
953 // parallel::shared::Triangulations, ownership of cells might change. If
954 // they allow artificial cells, we need to restore the "true" cell owners
955 // temporarily.
956 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
957 // the current set of subdomain ids, set subdomain ids to the "true" owner
958 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
959 // object, and later restore these flags when it is destroyed.
960 const ::internal::parallel::shared::
961 TemporarilyRestoreSubdomainIds<dim, spacedim>
962 subdomain_modifier(dof_handler->get_triangulation());
963
965 std::vector<types::global_dof_index> dofs;
966
967 typename std::map<std::pair<unsigned int, unsigned int>,
968 Pointerstruct>::const_iterator pointerstruct,
969 cell_map_end = cell_map.end();
970
971 Table<2, FullMatrix<double>> interpolation_hp;
972 internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
974
975 for (const auto &cell : dof_handler->cell_iterators())
976 {
977 pointerstruct =
978 cell_map.find(std::make_pair(cell->level(), cell->index()));
979
980 if (pointerstruct != cell_map_end)
981 {
982 const std::vector<types::global_dof_index> *const indexptr =
983 pointerstruct->second.indices_ptr;
984
985 const std::vector<Vector<typename VectorType::value_type>>
986 *const valuesptr = pointerstruct->second.dof_values_ptr;
987
988 // cell stayed as it was or was refined
989 if (indexptr != nullptr)
990 {
991 Assert(valuesptr == nullptr, ExcInternalError());
992
993 const unsigned int old_fe_index =
994 pointerstruct->second.active_fe_index;
995
996 // get the values of each of the input data vectors on this cell
997 // and prolong it to its children
998 unsigned int in_size = indexptr->size();
999 for (unsigned int j = 0; j < size; ++j)
1000 {
1001 tmp.reinit(in_size, true);
1002 for (unsigned int i = 0; i < in_size; ++i)
1004 all_in[j], (*indexptr)[i]);
1005
1006 cell->set_dof_values_by_interpolation(tmp,
1007 all_out[j],
1008 old_fe_index,
1009 true);
1010 }
1011 }
1012 else if (valuesptr)
1013 // the children of this cell were deleted
1014 {
1015 Assert(!cell->has_children(), ExcInternalError());
1016 Assert(indexptr == nullptr, ExcInternalError());
1017
1018 const unsigned int dofs_per_cell =
1019 cell->get_fe().n_dofs_per_cell();
1020 dofs.resize(dofs_per_cell);
1021 // get the local
1022 // indices
1023 cell->get_dof_indices(dofs);
1024
1025 // distribute the stored data to the new vectors
1026 for (unsigned int j = 0; j < size; ++j)
1027 {
1028 // make sure that the size of the stored indices is the same
1029 // as dofs_per_cell. this is kind of a test if we use the
1030 // same FE in the hp-case. to really do that test we would
1031 // have to store the fe_index of all cells
1033 nullptr;
1034 const unsigned int active_fe_index =
1035 cell->active_fe_index();
1036 if (active_fe_index !=
1037 pointerstruct->second.active_fe_index)
1038 {
1039 const unsigned int old_index =
1040 pointerstruct->second.active_fe_index;
1041 const FullMatrix<double> &interpolation_matrix =
1042 interpolation_hp(active_fe_index, old_index);
1043 // The interpolation matrix might be empty when using
1044 // FE_Nothing.
1045 if (interpolation_matrix.empty())
1046 tmp.reinit(dofs_per_cell, false);
1047 else
1048 {
1049 tmp.reinit(dofs_per_cell, true);
1050 AssertDimension((*valuesptr)[j].size(),
1051 interpolation_matrix.n());
1052 AssertDimension(tmp.size(),
1053 interpolation_matrix.m());
1054 interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
1055 }
1056 data = &tmp;
1057 }
1058 else
1059 data = &(*valuesptr)[j];
1060
1061
1062 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1064 (*data)(i), dofs[i], all_out[j]);
1065 }
1066 }
1067 // undefined status
1068 else
1069 Assert(false, ExcInternalError());
1070 }
1071 }
1072
1073 // We have written into the output vectors. If this was a PETSc vector, for
1074 // example, then we need to compress these to make future operations safe:
1075 for (auto &vec : all_out)
1076 vec.compress(VectorOperation::insert);
1077 }
1078
1079
1080
1081 template <int dim, typename VectorType, int spacedim>
1082 void
1084 const VectorType &in,
1085 VectorType &out) const
1086 {
1087 Assert(in.size() == n_dofs_old,
1088 ExcDimensionMismatch(in.size(), n_dofs_old));
1089 Assert(out.size() == dof_handler->n_dofs(),
1090 ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
1091
1092 std::vector<VectorType> all_in = {in};
1093 std::vector<VectorType> all_out = {out};
1094
1095 interpolate(all_in, all_out);
1096
1097 out = all_out[0];
1098 }
1099
1100
1101
1102 template <int dim, typename VectorType, int spacedim>
1103 std::size_t
1105 {
1106 // at the moment we do not include the memory
1107 // consumption of the cell_map as we have no
1108 // real idea about memory consumption of a
1109 // std::map
1110 return (MemoryConsumption::memory_consumption(dof_handler) +
1112 sizeof(prepared_for) +
1113 MemoryConsumption::memory_consumption(indices_on_cell) +
1114 MemoryConsumption::memory_consumption(dof_values_on_cell));
1115 }
1116
1117
1118
1119 template <int dim, typename VectorType, int spacedim>
1120 std::size_t
1122 memory_consumption() const
1123 {
1124 return sizeof(*this);
1125 }
1126
1127} // namespace Legacy
1128
1129
1130/*-------------- Explicit Instantiations -------------------------------*/
1131#define SPLIT_INSTANTIATIONS_COUNT 4
1132#ifndef SPLIT_INSTANTIATIONS_INDEX
1133# define SPLIT_INSTANTIATIONS_INDEX 0
1134#endif
1135#include "numerics/solution_transfer.inst"
1136
CellStatus
Definition cell_status.h:31
@ cell_will_be_refined
@ children_will_be_coarsened
const hp::FECollection< dim, spacedim > & get_fe_collection() const
bool has_hp_capabilities() const
size_type n() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
size_type m() const
void refine_interpolate(const VectorType &in, VectorType &out) const
void interpolate(const std::vector< VectorType > &all_in, std::vector< VectorType > &all_out) const
std::size_t memory_consumption() const
ObserverPointer< const DoFHandler< dim, spacedim >, SolutionTransfer< dim, VectorType, spacedim > > dof_handler
void prepare_for_coarsening_and_refinement(const std::vector< VectorType > &all_in)
SolutionTransfer(const DoFHandler< dim, spacedim > &dof)
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)
void deserialize(VectorType &in)
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)
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#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
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
void restriction_additive(const FiniteElement< dim, spacedim > &, std::vector< std::vector< bool > > &)
void extract_interpolation_matrices(const DoFHandler< dim, spacedim > &dof, ::Table< 2, FullMatrix< double > > &matrices)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)