Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2684-gc61376a70f 2025-02-22 15:30: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#ifdef DEBUG
372 for (const auto &child : cell->child_iterators())
373 Assert(child->is_active() && child->coarsen_flag_set(),
374 typename ::Triangulation<
375 dim>::ExcInconsistentCoarseningFlags());
376#endif
377
378 fe_index = ::internal::hp::DoFHandlerImplementation::
379 dominated_future_fe_on_children<dim, spacedim>(cell);
380 break;
381 }
382
383 default:
384 Assert(false, ExcInternalError());
385 break;
386 }
387 }
388
389 const unsigned int dofs_per_cell =
390 dof_handler->get_fe(fe_index).n_dofs_per_cell();
391
392 if (dofs_per_cell == 0)
393 return std::vector<char>(); // nothing to do for FE_Nothing
394
395 auto it_input = input_vectors.cbegin();
396 auto it_output = dof_values.begin();
397 for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
398 {
399 it_output->reinit(dofs_per_cell);
400 cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
401 }
402
403 return pack_dof_values<typename VectorType::value_type>(dof_values,
404 dofs_per_cell);
405}
406
407
408
409template <int dim, typename VectorType, int spacedim>
410void
413 const CellStatus status,
414 const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
415 std::vector<VectorType *> &all_out,
416 VectorType &valence)
417{
418 typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
419
420 unsigned int fe_index = 0;
421 if (dof_handler->has_hp_capabilities())
422 {
423 switch (status)
424 {
427 {
428 fe_index = cell->active_fe_index();
429 break;
430 }
431
433 {
434 // After refinement, this particular cell is no longer active,
435 // and its children have inherited its FE index. However, to
436 // unpack the data on the old cell, we need to recover its FE
437 // index from one of the children. Just to be sure, we also
438 // check if all children have the same FE index.
439 fe_index = cell->child(0)->active_fe_index();
440 for (unsigned int child_index = 1;
441 child_index < cell->n_children();
442 ++child_index)
443 Assert(cell->child(child_index)->active_fe_index() == fe_index,
445 break;
446 }
447
448 default:
449 Assert(false, ExcInternalError());
450 break;
451 }
452 }
453
454 const unsigned int dofs_per_cell =
455 dof_handler->get_fe(fe_index).n_dofs_per_cell();
456
457 if (dofs_per_cell == 0)
458 return; // nothing to do for FE_Nothing
459
460 const std::vector<::Vector<typename VectorType::value_type>>
461 dof_values =
462 unpack_dof_values<typename VectorType::value_type>(data_range,
463 dofs_per_cell);
464
465 // check if sizes match
466 AssertDimension(dof_values.size(), all_out.size());
467
468 // check if we have enough dofs provided by the FE object
469 // to interpolate the transferred data correctly
470 for (auto it_dof_values = dof_values.begin();
471 it_dof_values != dof_values.end();
472 ++it_dof_values)
473 Assert(
474 dofs_per_cell == it_dof_values->size(),
476 "The transferred data was packed with a different number of dofs than the "
477 "currently registered FE object assigned to the DoFHandler has."));
478
479 // distribute data for each registered vector on mesh
480 auto it_input = dof_values.cbegin();
481 auto it_output = all_out.begin();
482 for (; it_input != dof_values.cend(); ++it_input, ++it_output)
483 if (average_values)
484 cell->distribute_local_to_global_by_interpolation(*it_input,
485 *(*it_output),
486 fe_index);
487 else
488 cell->set_dof_values_by_interpolation(*it_input,
489 *(*it_output),
490 fe_index,
491 true);
492
493 if (average_values)
494 {
495 // compute valence vector if averaging should be performed
496 Vector<typename VectorType::value_type> ones(dofs_per_cell);
497 ones = 1.0;
498 cell->distribute_local_to_global_by_interpolation(ones,
499 valence,
500 fe_index);
501 }
502}
503
504
505
506template <int dim, typename VectorType, int spacedim>
507void
512
513
514
515namespace Legacy
516{
517
518 template <int dim, typename VectorType, int spacedim>
520 const DoFHandler<dim, spacedim> &dof)
521 : dof_handler(&dof, typeid(*this).name())
522 , n_dofs_old(0)
523 , prepared_for(none)
524 {
525 Assert((dynamic_cast<
527 &dof_handler->get_triangulation()) == nullptr),
529 "You are calling the ::SolutionTransfer class "
530 "with a DoFHandler that is built on a "
531 "parallel::distributed::Triangulation. This will not "
532 "work for parallel computations. You probably want to "
533 "use the parallel::distributed::SolutionTransfer class."));
534 }
535
536
537
538 template <int dim, typename VectorType, int spacedim>
543
544
545
546 template <int dim, typename VectorType, int spacedim>
547 void
549 {
550 indices_on_cell.clear();
551 dof_values_on_cell.clear();
552 cell_map.clear();
553
554 prepared_for = none;
555 }
556
557
558
559 template <int dim, typename VectorType, int spacedim>
560 void
562 {
563 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
564 Assert(prepared_for != coarsening_and_refinement,
565 ExcAlreadyPrepForCoarseAndRef());
566
567 clear();
568
569 // We need to access dof indices on the entire domain. For
570 // parallel::shared::Triangulations, ownership of cells might change. If
571 // they allow artificial cells, we need to restore the "true" cell owners
572 // temporarily.
573 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
574 // the current set of subdomain ids, set subdomain ids to the "true" owner
575 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
576 // object, and later restore these flags when it is destroyed.
577 const ::internal::parallel::shared::
578 TemporarilyRestoreSubdomainIds<dim, spacedim>
579 subdomain_modifier(dof_handler->get_triangulation());
580
581 const unsigned int n_active_cells =
582 dof_handler->get_triangulation().n_active_cells();
583 n_dofs_old = dof_handler->n_dofs();
584
585 // efficient reallocation of indices_on_cell
586 std::vector<std::vector<types::global_dof_index>>(n_active_cells)
587 .swap(indices_on_cell);
588
589 for (const auto &cell : dof_handler->active_cell_iterators())
590 {
591 const unsigned int i = cell->active_cell_index();
592 indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
593 // on each cell store the indices of the
594 // dofs. after refining we get the values
595 // on the children by taking these
596 // indices, getting the respective values
597 // out of the data vectors and prolonging
598 // them to the children
599 cell->get_dof_indices(indices_on_cell[i]);
600 cell_map[std::make_pair(cell->level(), cell->index())] =
601 Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
602 }
603 prepared_for = pure_refinement;
604 }
605
606
607
608 template <int dim, typename VectorType, int spacedim>
609 void
611 const VectorType &in,
612 VectorType &out) const
613 {
614 Assert(prepared_for == pure_refinement, ExcNotPrepared());
615 Assert(in.size() == n_dofs_old,
616 ExcDimensionMismatch(in.size(), n_dofs_old));
617 Assert(out.size() == dof_handler->n_dofs(),
618 ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
619 Assert(&in != &out,
620 ExcMessage("Vectors cannot be used as input and output"
621 " at the same time!"));
622
623 // We need to access dof indices on the entire domain. For
624 // parallel::shared::Triangulations, ownership of cells might change. If
625 // they allow artificial cells, we need to restore the "true" cell owners
626 // temporarily.
627 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
628 // the current set of subdomain ids, set subdomain ids to the "true" owner
629 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
630 // object, and later restore these flags when it is destroyed.
631 const ::internal::parallel::shared::
632 TemporarilyRestoreSubdomainIds<dim, spacedim>
633 subdomain_modifier(dof_handler->get_triangulation());
634
636
637 typename std::map<std::pair<unsigned int, unsigned int>,
638 Pointerstruct>::const_iterator pointerstruct,
639 cell_map_end = cell_map.end();
640
641 for (const auto &cell : dof_handler->cell_iterators())
642 {
643 pointerstruct =
644 cell_map.find(std::make_pair(cell->level(), cell->index()));
645
646 if (pointerstruct != cell_map_end)
647 // this cell was refined or not
648 // touched at all, so we can get
649 // the new values by just setting
650 // or interpolating to the children,
651 // which is both done by one
652 // function
653 {
654 const unsigned int this_fe_index =
655 pointerstruct->second.active_fe_index;
656 const unsigned int dofs_per_cell =
657 cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
658 local_values.reinit(dofs_per_cell, true);
659
660 // make sure that the size of the stored indices is the same as
661 // dofs_per_cell. since we store the desired fe_index, we know
662 // what this size should be
663 Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
665 for (unsigned int i = 0; i < dofs_per_cell; ++i)
666 local_values(i) =
668 in, (*pointerstruct->second.indices_ptr)[i]);
669 cell->set_dof_values_by_interpolation(local_values,
670 out,
671 this_fe_index,
672 true);
673 }
674 }
675 }
676
677
678
679 namespace internal
680 {
694 template <int dim, int spacedim>
695 void
697 const DoFHandler<dim, spacedim> &dof,
698 ::Table<2, FullMatrix<double>> &matrices)
699 {
700 if (dof.has_hp_capabilities() == false)
701 return;
702
703 const ::hp::FECollection<dim, spacedim> &fe =
704 dof.get_fe_collection();
705 matrices.reinit(fe.size(), fe.size());
706 for (unsigned int i = 0; i < fe.size(); ++i)
707 for (unsigned int j = 0; j < fe.size(); ++j)
708 if (i != j)
709 {
710 matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
711 fe[j].n_dofs_per_cell());
712
713 // see if we can get the interpolation matrices for this
714 // combination of elements. if not, reset the matrix sizes to zero
715 // to indicate that this particular combination isn't
716 // supported. this isn't an outright error right away since we may
717 // never need to actually interpolate between these two elements
718 // on actual cells; we simply have to trigger an error if someone
719 // actually tries
720 try
721 {
722 fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
723 }
724 catch (const typename FiniteElement<dim, spacedim>::
725 ExcInterpolationNotImplemented &)
726 {
727 matrices(i, j).reinit(0, 0);
728 }
729 }
730 }
731
732
733 template <int dim, int spacedim>
734 void
736 std::vector<std::vector<bool>> &)
737 {}
738
739 template <int dim, int spacedim>
740 void
742 const ::hp::FECollection<dim, spacedim> &fe,
743 std::vector<std::vector<bool>> &restriction_is_additive)
744 {
745 restriction_is_additive.resize(fe.size());
746 for (unsigned int f = 0; f < fe.size(); ++f)
747 {
748 restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
749 for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
750 restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
751 }
752 }
753 } // namespace internal
754
755
756
757 template <int dim, typename VectorType, int spacedim>
758 void
760 prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
761 {
762 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
763 Assert(prepared_for != coarsening_and_refinement,
764 ExcAlreadyPrepForCoarseAndRef());
765
766 clear();
767 n_dofs_old = dof_handler->n_dofs();
768 const unsigned int in_size = all_in.size();
769
770#ifdef DEBUG
771 Assert(in_size != 0,
772 ExcMessage("The array of input vectors you pass to this "
773 "function has no elements. This is not useful."));
774 for (unsigned int i = 0; i < in_size; ++i)
775 {
776 Assert(all_in[i].size() == n_dofs_old,
777 ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
778 }
779#endif
780
781 // We need to access dof indices on the entire domain. For
782 // parallel::shared::Triangulations, ownership of cells might change. If
783 // they allow artificial cells, we need to restore the "true" cell owners
784 // temporarily.
785 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
786 // the current set of subdomain ids, set subdomain ids to the "true" owner
787 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
788 // object, and later restore these flags when it is destroyed.
789 const ::internal::parallel::shared::
790 TemporarilyRestoreSubdomainIds<dim, spacedim>
791 subdomain_modifier(dof_handler->get_triangulation());
792
793 // first count the number
794 // of cells that will be coarsened
795 // and that'll stay or be refined
796 unsigned int n_cells_to_coarsen = 0;
797 unsigned int n_cells_to_stay_or_refine = 0;
798 for (const auto &act_cell : dof_handler->active_cell_iterators())
799 {
800 if (act_cell->coarsen_flag_set())
801 ++n_cells_to_coarsen;
802 else
803 ++n_cells_to_stay_or_refine;
804 }
805 Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
806 dof_handler->get_triangulation().n_active_cells(),
808
809 unsigned int n_coarsen_fathers = 0;
810 for (const auto &cell : dof_handler->cell_iterators())
811 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
812 ++n_coarsen_fathers;
813 Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
814 (void)n_cells_to_coarsen;
815
816 // allocate the needed memory. initialize
817 // the following arrays in an efficient
818 // way, without copying much
819 std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
820 .swap(indices_on_cell);
821
822 std::vector<std::vector<Vector<typename VectorType::value_type>>>(
823 n_coarsen_fathers,
824 std::vector<Vector<typename VectorType::value_type>>(in_size))
825 .swap(dof_values_on_cell);
826
827 Table<2, FullMatrix<double>> interpolation_hp;
828 std::vector<std::vector<bool>> restriction_is_additive;
829
830 internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
831 internal::restriction_additive(dof_handler->get_fe_collection(),
832 restriction_is_additive);
833
834 // we need counters for
835 // the 'to_stay_or_refine' cells 'n_sr' and
836 // the 'coarsen_fathers' cells 'n_cf',
837 unsigned int n_sr = 0, n_cf = 0;
838 for (const auto &cell : dof_handler->cell_iterators())
839 {
840 // CASE 1: active cell that remains as it is
841 if (cell->is_active() && !cell->coarsen_flag_set())
842 {
843 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
844 indices_on_cell[n_sr].resize(dofs_per_cell);
845 // cell will not be coarsened,
846 // so we get away by storing the
847 // dof indices and later
848 // interpolating to the children
849 cell->get_dof_indices(indices_on_cell[n_sr]);
850 cell_map[std::make_pair(cell->level(), cell->index())] =
851 Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
852 ++n_sr;
853 }
854
855 // CASE 2: cell is inactive but will become active
856 else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
857 {
858 // we will need to interpolate from the children of this cell
859 // to the current one. in the hp-context, this also means
860 // we need to figure out which finite element space to interpolate
861 // to since that is not implied by the global FE as in the non-hp-
862 // case. we choose the 'least dominant fe' on all children from
863 // the associated FECollection.
864 std::set<unsigned int> fe_indices_children;
865 for (const auto &child : cell->child_iterators())
866 {
867 Assert(child->is_active() && child->coarsen_flag_set(),
868 typename ::Triangulation<
869 dim>::ExcInconsistentCoarseningFlags());
870
871 fe_indices_children.insert(child->active_fe_index());
872 }
873 Assert(!fe_indices_children.empty(), ExcInternalError());
874
875 const unsigned int target_fe_index =
876 dof_handler->get_fe_collection().find_dominated_fe_extended(
877 fe_indices_children, /*codim=*/0);
878
879 Assert(target_fe_index != numbers::invalid_unsigned_int,
880 ::internal::hp::DoFHandlerImplementation::
881 ExcNoDominatedFiniteElementOnChildren());
882
883 const unsigned int dofs_per_cell =
884 dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
885
886 std::vector<Vector<typename VectorType::value_type>>(
887 in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
888 .swap(dof_values_on_cell[n_cf]);
889
890
891 // store the data of each of the input vectors. get this data
892 // as interpolated onto a finite element space that encompasses
893 // that of all the children. note that
894 // cell->get_interpolated_dof_values already does all of the
895 // interpolations between spaces
896 for (unsigned int j = 0; j < in_size; ++j)
897 cell->get_interpolated_dof_values(all_in[j],
898 dof_values_on_cell[n_cf][j],
899 target_fe_index);
900 cell_map[std::make_pair(cell->level(), cell->index())] =
901 Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
902 ++n_cf;
903 }
904 }
905 Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
906 Assert(n_cf == n_coarsen_fathers, ExcInternalError());
907
908 prepared_for = coarsening_and_refinement;
909 }
910
911
912
913 template <int dim, typename VectorType, int spacedim>
914 void
916 prepare_for_coarsening_and_refinement(const VectorType &in)
917 {
918 std::vector<VectorType> all_in(1, in);
919 prepare_for_coarsening_and_refinement(all_in);
920 }
921
922
923
924 template <int dim, typename VectorType, int spacedim>
925 void
927 const std::vector<VectorType> &all_in,
928 std::vector<VectorType> &all_out) const
929 {
930 const unsigned int size = all_in.size();
931#ifdef DEBUG
932 Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
933 Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
934 for (unsigned int i = 0; i < size; ++i)
935 Assert(all_in[i].size() == n_dofs_old,
936 ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
937 for (unsigned int i = 0; i < all_out.size(); ++i)
938 Assert(all_out[i].size() == dof_handler->n_dofs(),
939 ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
940 for (unsigned int i = 0; i < size; ++i)
941 for (unsigned int j = 0; j < size; ++j)
942 Assert(&all_in[i] != &all_out[j],
943 ExcMessage("Vectors cannot be used as input and output"
944 " at the same time!"));
945#endif
946
947 // We need to access dof indices on the entire domain. For
948 // parallel::shared::Triangulations, ownership of cells might change. If
949 // they allow artificial cells, we need to restore the "true" cell owners
950 // temporarily.
951 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
952 // the current set of subdomain ids, set subdomain ids to the "true" owner
953 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
954 // object, and later restore these flags when it is destroyed.
955 const ::internal::parallel::shared::
956 TemporarilyRestoreSubdomainIds<dim, spacedim>
957 subdomain_modifier(dof_handler->get_triangulation());
958
960 std::vector<types::global_dof_index> dofs;
961
962 typename std::map<std::pair<unsigned int, unsigned int>,
963 Pointerstruct>::const_iterator pointerstruct,
964 cell_map_end = cell_map.end();
965
966 Table<2, FullMatrix<double>> interpolation_hp;
967 internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
969
970 for (const auto &cell : dof_handler->cell_iterators())
971 {
972 pointerstruct =
973 cell_map.find(std::make_pair(cell->level(), cell->index()));
974
975 if (pointerstruct != cell_map_end)
976 {
977 const std::vector<types::global_dof_index> *const indexptr =
978 pointerstruct->second.indices_ptr;
979
980 const std::vector<Vector<typename VectorType::value_type>>
981 *const valuesptr = pointerstruct->second.dof_values_ptr;
982
983 // cell stayed as it was or was refined
984 if (indexptr != nullptr)
985 {
986 Assert(valuesptr == nullptr, ExcInternalError());
987
988 const unsigned int old_fe_index =
989 pointerstruct->second.active_fe_index;
990
991 // get the values of each of the input data vectors on this cell
992 // and prolong it to its children
993 unsigned int in_size = indexptr->size();
994 for (unsigned int j = 0; j < size; ++j)
995 {
996 tmp.reinit(in_size, true);
997 for (unsigned int i = 0; i < in_size; ++i)
999 all_in[j], (*indexptr)[i]);
1000
1001 cell->set_dof_values_by_interpolation(tmp,
1002 all_out[j],
1003 old_fe_index,
1004 true);
1005 }
1006 }
1007 else if (valuesptr)
1008 // the children of this cell were deleted
1009 {
1010 Assert(!cell->has_children(), ExcInternalError());
1011 Assert(indexptr == nullptr, ExcInternalError());
1012
1013 const unsigned int dofs_per_cell =
1014 cell->get_fe().n_dofs_per_cell();
1015 dofs.resize(dofs_per_cell);
1016 // get the local
1017 // indices
1018 cell->get_dof_indices(dofs);
1019
1020 // distribute the stored data to the new vectors
1021 for (unsigned int j = 0; j < size; ++j)
1022 {
1023 // make sure that the size of the stored indices is the same
1024 // as dofs_per_cell. this is kind of a test if we use the
1025 // same FE in the hp-case. to really do that test we would
1026 // have to store the fe_index of all cells
1028 nullptr;
1029 const unsigned int active_fe_index =
1030 cell->active_fe_index();
1031 if (active_fe_index !=
1032 pointerstruct->second.active_fe_index)
1033 {
1034 const unsigned int old_index =
1035 pointerstruct->second.active_fe_index;
1036 const FullMatrix<double> &interpolation_matrix =
1037 interpolation_hp(active_fe_index, old_index);
1038 // The interpolation matrix might be empty when using
1039 // FE_Nothing.
1040 if (interpolation_matrix.empty())
1041 tmp.reinit(dofs_per_cell, false);
1042 else
1043 {
1044 tmp.reinit(dofs_per_cell, true);
1045 AssertDimension((*valuesptr)[j].size(),
1046 interpolation_matrix.n());
1047 AssertDimension(tmp.size(),
1048 interpolation_matrix.m());
1049 interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
1050 }
1051 data = &tmp;
1052 }
1053 else
1054 data = &(*valuesptr)[j];
1055
1056
1057 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1059 (*data)(i), dofs[i], all_out[j]);
1060 }
1061 }
1062 // undefined status
1063 else
1064 Assert(false, ExcInternalError());
1065 }
1066 }
1067
1068 // We have written into the output vectors. If this was a PETSc vector, for
1069 // example, then we need to compress these to make future operations safe:
1070 for (auto &vec : all_out)
1071 vec.compress(VectorOperation::insert);
1072 }
1073
1074
1075
1076 template <int dim, typename VectorType, int spacedim>
1077 void
1079 const VectorType &in,
1080 VectorType &out) const
1081 {
1082 Assert(in.size() == n_dofs_old,
1083 ExcDimensionMismatch(in.size(), n_dofs_old));
1084 Assert(out.size() == dof_handler->n_dofs(),
1085 ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
1086
1087 std::vector<VectorType> all_in = {in};
1088 std::vector<VectorType> all_out = {out};
1089
1090 interpolate(all_in, all_out);
1091
1092 out = all_out[0];
1093 }
1094
1095
1096
1097 template <int dim, typename VectorType, int spacedim>
1098 std::size_t
1100 {
1101 // at the moment we do not include the memory
1102 // consumption of the cell_map as we have no
1103 // real idea about memory consumption of a
1104 // std::map
1105 return (MemoryConsumption::memory_consumption(dof_handler) +
1107 sizeof(prepared_for) +
1108 MemoryConsumption::memory_consumption(indices_on_cell) +
1109 MemoryConsumption::memory_consumption(dof_values_on_cell));
1110 }
1111
1112
1113
1114 template <int dim, typename VectorType, int spacedim>
1115 std::size_t
1117 memory_consumption() const
1118 {
1119 return sizeof(*this);
1120 }
1121
1122} // namespace Legacy
1123
1124
1125/*-------------- Explicit Instantiations -------------------------------*/
1126#define SPLIT_INSTANTIATIONS_COUNT 4
1127#ifndef SPLIT_INSTANTIATIONS_INDEX
1128# define SPLIT_INSTANTIATIONS_INDEX 0
1129#endif
1130#include "numerics/solution_transfer.inst"
1131
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:522
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:523
#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:740
std::size_t size
Definition mpi.cc:739
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)