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