Reference documentation for deal.II version GIT relicensing-43-g3fc2fae7c3 2024-03-03 21:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
solution_transfer.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 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
16
19
22
23#include <deal.II/fe/fe.h>
24
25#include <deal.II/grid/tria.h>
28
36#include <deal.II/lac/vector.h>
38
40
42
43template <int dim, typename VectorType, int spacedim>
46 : dof_handler(&dof, typeid(*this).name())
47 , n_dofs_old(0)
48 , prepared_for(none)
49{
50 Assert(
52 &dof_handler->get_triangulation()) == nullptr),
53 ExcMessage("You are calling the ::SolutionTransfer class "
54 "with a DoFHandler that is built on a "
55 "parallel::distributed::Triangulation. This will not "
56 "work for parallel computations. You probably want to "
57 "use the parallel::distributed::SolutionTransfer class."));
58}
59
60
61
62template <int dim, typename VectorType, int spacedim>
67
68
69
70template <int dim, typename VectorType, int spacedim>
71void
73{
74 indices_on_cell.clear();
75 dof_values_on_cell.clear();
76 cell_map.clear();
77
78 prepared_for = none;
79}
80
81
82
83template <int dim, typename VectorType, int spacedim>
84void
86{
87 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
88 Assert(prepared_for != coarsening_and_refinement,
89 ExcAlreadyPrepForCoarseAndRef());
90
91 clear();
92
93 // We need to access dof indices on the entire domain. For
94 // parallel::shared::Triangulations, ownership of cells might change. If they
95 // allow artificial cells, we need to restore the "true" cell owners
96 // temporarily.
97 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
98 // the current set of subdomain ids, set subdomain ids to the "true" owner of
99 // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
100 // and later restore these flags when it is destroyed.
102 spacedim>
103 subdomain_modifier(dof_handler->get_triangulation());
104
105 const unsigned int n_active_cells =
106 dof_handler->get_triangulation().n_active_cells();
107 n_dofs_old = dof_handler->n_dofs();
108
109 // efficient reallocation of indices_on_cell
110 std::vector<std::vector<types::global_dof_index>>(n_active_cells)
111 .swap(indices_on_cell);
112
113 for (const auto &cell : dof_handler->active_cell_iterators())
114 {
115 const unsigned int i = cell->active_cell_index();
116 indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
117 // on each cell store the indices of the
118 // dofs. after refining we get the values
119 // on the children by taking these
120 // indices, getting the respective values
121 // out of the data vectors and prolonging
122 // them to the children
123 cell->get_dof_indices(indices_on_cell[i]);
124 cell_map[std::make_pair(cell->level(), cell->index())] =
125 Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
126 }
127 prepared_for = pure_refinement;
128}
129
130
131
132template <int dim, typename VectorType, int spacedim>
133void
135 const VectorType &in,
136 VectorType &out) const
137{
138 Assert(prepared_for == pure_refinement, ExcNotPrepared());
139 Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
140 Assert(out.size() == dof_handler->n_dofs(),
141 ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
142 Assert(&in != &out,
143 ExcMessage("Vectors cannot be used as input and output"
144 " at the same time!"));
145
146 // We need to access dof indices on the entire domain. For
147 // parallel::shared::Triangulations, ownership of cells might change. If they
148 // allow artificial cells, we need to restore the "true" cell owners
149 // temporarily.
150 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
151 // the current set of subdomain ids, set subdomain ids to the "true" owner of
152 // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
153 // and later restore these flags when it is destroyed.
155 spacedim>
156 subdomain_modifier(dof_handler->get_triangulation());
157
159
160 typename std::map<std::pair<unsigned int, unsigned int>,
161 Pointerstruct>::const_iterator pointerstruct,
162 cell_map_end = cell_map.end();
163
164 for (const auto &cell : dof_handler->cell_iterators())
165 {
166 pointerstruct =
167 cell_map.find(std::make_pair(cell->level(), cell->index()));
168
169 if (pointerstruct != cell_map_end)
170 // this cell was refined or not
171 // touched at all, so we can get
172 // the new values by just setting
173 // or interpolating to the children,
174 // which is both done by one
175 // function
176 {
177 const unsigned int this_fe_index =
178 pointerstruct->second.active_fe_index;
179 const unsigned int dofs_per_cell =
180 cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
181 local_values.reinit(dofs_per_cell, true);
182
183 // make sure that the size of the stored indices is the same as
184 // dofs_per_cell. since we store the desired fe_index, we know
185 // what this size should be
186 Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
188 for (unsigned int i = 0; i < dofs_per_cell; ++i)
190 in, (*pointerstruct->second.indices_ptr)[i]);
191 cell->set_dof_values_by_interpolation(local_values,
192 out,
193 this_fe_index,
194 true);
195 }
196 }
197}
198
199
200
201namespace internal
202{
216 template <int dim, int spacedim>
217 void
219 ::Table<2, FullMatrix<double>> &matrices)
220 {
221 if (dof.has_hp_capabilities() == false)
222 return;
223
224 const ::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
225 matrices.reinit(fe.size(), fe.size());
226 for (unsigned int i = 0; i < fe.size(); ++i)
227 for (unsigned int j = 0; j < fe.size(); ++j)
228 if (i != j)
229 {
230 matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
231 fe[j].n_dofs_per_cell());
232
233 // see if we can get the interpolation matrices for this
234 // combination of elements. if not, reset the matrix sizes to zero
235 // to indicate that this particular combination isn't
236 // supported. this isn't an outright error right away since we may
237 // never need to actually interpolate between these two elements
238 // on actual cells; we simply have to trigger an error if someone
239 // actually tries
240 try
241 {
242 fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
243 }
244 catch (const typename FiniteElement<dim, spacedim>::
245 ExcInterpolationNotImplemented &)
246 {
247 matrices(i, j).reinit(0, 0);
248 }
249 }
250 }
251
252
253 template <int dim, int spacedim>
254 void
256 std::vector<std::vector<bool>> &)
257 {}
258
259 template <int dim, int spacedim>
260 void
261 restriction_additive(const ::hp::FECollection<dim, spacedim> &fe,
262 std::vector<std::vector<bool>> &restriction_is_additive)
263 {
264 restriction_is_additive.resize(fe.size());
265 for (unsigned int f = 0; f < fe.size(); ++f)
266 {
267 restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
268 for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
269 restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
270 }
271 }
272} // namespace internal
273
274
275
276template <int dim, typename VectorType, int spacedim>
277void
279 prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
280{
281 Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
282 Assert(prepared_for != coarsening_and_refinement,
283 ExcAlreadyPrepForCoarseAndRef());
284
285 clear();
286 n_dofs_old = dof_handler->n_dofs();
287 const unsigned int in_size = all_in.size();
288
289#ifdef DEBUG
290 Assert(in_size != 0,
291 ExcMessage("The array of input vectors you pass to this "
292 "function has no elements. This is not useful."));
293 for (unsigned int i = 0; i < in_size; ++i)
294 {
295 Assert(all_in[i].size() == n_dofs_old,
296 ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
297 }
298#endif
299
300 // We need to access dof indices on the entire domain. For
301 // parallel::shared::Triangulations, ownership of cells might change. If they
302 // allow artificial cells, we need to restore the "true" cell owners
303 // temporarily.
304 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
305 // the current set of subdomain ids, set subdomain ids to the "true" owner of
306 // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
307 // and later restore these flags when it is destroyed.
309 spacedim>
310 subdomain_modifier(dof_handler->get_triangulation());
311
312 // first count the number
313 // of cells that will be coarsened
314 // and that'll stay or be refined
315 unsigned int n_cells_to_coarsen = 0;
316 unsigned int n_cells_to_stay_or_refine = 0;
317 for (const auto &act_cell : dof_handler->active_cell_iterators())
318 {
319 if (act_cell->coarsen_flag_set())
320 ++n_cells_to_coarsen;
321 else
322 ++n_cells_to_stay_or_refine;
323 }
324 Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
325 dof_handler->get_triangulation().n_active_cells(),
327
328 unsigned int n_coarsen_fathers = 0;
329 for (const auto &cell : dof_handler->cell_iterators())
330 if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
331 ++n_coarsen_fathers;
332 Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
333 (void)n_cells_to_coarsen;
334
335 // allocate the needed memory. initialize
336 // the following arrays in an efficient
337 // way, without copying much
338 std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
339 .swap(indices_on_cell);
340
341 std::vector<std::vector<Vector<typename VectorType::value_type>>>(
342 n_coarsen_fathers,
343 std::vector<Vector<typename VectorType::value_type>>(in_size))
344 .swap(dof_values_on_cell);
345
346 Table<2, FullMatrix<double>> interpolation_hp;
347 std::vector<std::vector<bool>> restriction_is_additive;
348
349 internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
350 internal::restriction_additive(dof_handler->get_fe_collection(),
351 restriction_is_additive);
352
353 // we need counters for
354 // the 'to_stay_or_refine' cells 'n_sr' and
355 // the 'coarsen_fathers' cells 'n_cf',
356 unsigned int n_sr = 0, n_cf = 0;
357 for (const auto &cell : dof_handler->cell_iterators())
358 {
359 // CASE 1: active cell that remains as it is
360 if (cell->is_active() && !cell->coarsen_flag_set())
361 {
362 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
363 indices_on_cell[n_sr].resize(dofs_per_cell);
364 // cell will not be coarsened,
365 // so we get away by storing the
366 // dof indices and later
367 // interpolating to the children
368 cell->get_dof_indices(indices_on_cell[n_sr]);
369 cell_map[std::make_pair(cell->level(), cell->index())] =
370 Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
371 ++n_sr;
372 }
373
374 // CASE 2: cell is inactive but will become active
375 else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
376 {
377 // we will need to interpolate from the children of this cell
378 // to the current one. in the hp-context, this also means
379 // we need to figure out which finite element space to interpolate
380 // to since that is not implied by the global FE as in the non-hp-
381 // case. we choose the 'least dominant fe' on all children from
382 // the associated FECollection.
383 std::set<unsigned int> fe_indices_children;
384 for (const auto &child : cell->child_iterators())
385 {
386 Assert(child->is_active() && child->coarsen_flag_set(),
387 typename ::Triangulation<
388 dim>::ExcInconsistentCoarseningFlags());
389
390 fe_indices_children.insert(child->active_fe_index());
391 }
392 Assert(!fe_indices_children.empty(), ExcInternalError());
393
394 const unsigned int target_fe_index =
395 dof_handler->get_fe_collection().find_dominated_fe_extended(
396 fe_indices_children, /*codim=*/0);
397
398 Assert(target_fe_index != numbers::invalid_unsigned_int,
400 ExcNoDominatedFiniteElementOnChildren());
401
402 const unsigned int dofs_per_cell =
403 dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
404
405 std::vector<Vector<typename VectorType::value_type>>(
406 in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
407 .swap(dof_values_on_cell[n_cf]);
408
409
410 // store the data of each of the input vectors. get this data
411 // as interpolated onto a finite element space that encompasses
412 // that of all the children. note that
413 // cell->get_interpolated_dof_values already does all of the
414 // interpolations between spaces
415 for (unsigned int j = 0; j < in_size; ++j)
416 cell->get_interpolated_dof_values(all_in[j],
417 dof_values_on_cell[n_cf][j],
418 target_fe_index);
419 cell_map[std::make_pair(cell->level(), cell->index())] =
420 Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
421 ++n_cf;
422 }
423 }
424 Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
425 Assert(n_cf == n_coarsen_fathers, ExcInternalError());
426
427 prepared_for = coarsening_and_refinement;
428}
429
430
431
432template <int dim, typename VectorType, int spacedim>
433void
435 prepare_for_coarsening_and_refinement(const VectorType &in)
436{
437 std::vector<VectorType> all_in(1, in);
438 prepare_for_coarsening_and_refinement(all_in);
439}
440
441
442
443template <int dim, typename VectorType, int spacedim>
444void
446 const std::vector<VectorType> &all_in,
447 std::vector<VectorType> &all_out) const
448{
449 const unsigned int size = all_in.size();
450#ifdef DEBUG
451 Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
452 Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
453 for (unsigned int i = 0; i < size; ++i)
454 Assert(all_in[i].size() == n_dofs_old,
455 ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
456 for (unsigned int i = 0; i < all_out.size(); ++i)
457 Assert(all_out[i].size() == dof_handler->n_dofs(),
458 ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
459 for (unsigned int i = 0; i < size; ++i)
460 for (unsigned int j = 0; j < size; ++j)
461 Assert(&all_in[i] != &all_out[j],
462 ExcMessage("Vectors cannot be used as input and output"
463 " at the same time!"));
464#endif
465
466 // We need to access dof indices on the entire domain. For
467 // parallel::shared::Triangulations, ownership of cells might change. If they
468 // allow artificial cells, we need to restore the "true" cell owners
469 // temporarily.
470 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
471 // the current set of subdomain ids, set subdomain ids to the "true" owner of
472 // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
473 // and later restore these flags when it is destroyed.
475 spacedim>
476 subdomain_modifier(dof_handler->get_triangulation());
477
479 std::vector<types::global_dof_index> dofs;
480
481 typename std::map<std::pair<unsigned int, unsigned int>,
482 Pointerstruct>::const_iterator pointerstruct,
483 cell_map_end = cell_map.end();
484
485 Table<2, FullMatrix<double>> interpolation_hp;
486 internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
488
489 for (const auto &cell : dof_handler->cell_iterators())
490 {
491 pointerstruct =
492 cell_map.find(std::make_pair(cell->level(), cell->index()));
493
494 if (pointerstruct != cell_map_end)
495 {
496 const std::vector<types::global_dof_index> *const indexptr =
497 pointerstruct->second.indices_ptr;
498
499 const std::vector<Vector<typename VectorType::value_type>>
500 *const valuesptr = pointerstruct->second.dof_values_ptr;
501
502 // cell stayed as it was or was refined
503 if (indexptr != nullptr)
504 {
505 Assert(valuesptr == nullptr, ExcInternalError());
506
507 const unsigned int old_fe_index =
508 pointerstruct->second.active_fe_index;
509
510 // get the values of each of the input data vectors on this cell
511 // and prolong it to its children
512 unsigned int in_size = indexptr->size();
513 for (unsigned int j = 0; j < size; ++j)
514 {
515 tmp.reinit(in_size, true);
516 for (unsigned int i = 0; i < in_size; ++i)
517 tmp(i) =
519 (*indexptr)[i]);
520
521 cell->set_dof_values_by_interpolation(tmp,
522 all_out[j],
523 old_fe_index,
524 true);
525 }
526 }
527 else if (valuesptr)
528 // the children of this cell were deleted
529 {
530 Assert(!cell->has_children(), ExcInternalError());
531 Assert(indexptr == nullptr, ExcInternalError());
532
533 const unsigned int dofs_per_cell =
534 cell->get_fe().n_dofs_per_cell();
535 dofs.resize(dofs_per_cell);
536 // get the local
537 // indices
538 cell->get_dof_indices(dofs);
539
540 // distribute the stored data to the new vectors
541 for (unsigned int j = 0; j < size; ++j)
542 {
543 // make sure that the size of the stored indices is the same
544 // as dofs_per_cell. this is kind of a test if we use the same
545 // FE in the hp-case. to really do that test we would have to
546 // store the fe_index of all cells
547 const Vector<typename VectorType::value_type> *data = nullptr;
548 const unsigned int active_fe_index = cell->active_fe_index();
549 if (active_fe_index != pointerstruct->second.active_fe_index)
550 {
551 const unsigned int old_index =
552 pointerstruct->second.active_fe_index;
553 const FullMatrix<double> &interpolation_matrix =
554 interpolation_hp(active_fe_index, old_index);
555 // The interpolation matrix might be empty when using
556 // FE_Nothing.
557 if (interpolation_matrix.empty())
558 tmp.reinit(dofs_per_cell, false);
559 else
560 {
561 tmp.reinit(dofs_per_cell, true);
562 AssertDimension((*valuesptr)[j].size(),
563 interpolation_matrix.n());
564 AssertDimension(tmp.size(), interpolation_matrix.m());
565 interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
566 }
567 data = &tmp;
568 }
569 else
570 data = &(*valuesptr)[j];
571
572
573 for (unsigned int i = 0; i < dofs_per_cell; ++i)
575 dofs[i],
576 all_out[j]);
577 }
578 }
579 // undefined status
580 else
582 }
583 }
584
585 // We have written into the output vectors. If this was a PETSc vector, for
586 // example, then we need to compress these to make future operations safe:
587 for (auto &vec : all_out)
588 vec.compress(VectorOperation::insert);
589}
590
591
592
593template <int dim, typename VectorType, int spacedim>
594void
596 VectorType &out) const
597{
598 Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
599 Assert(out.size() == dof_handler->n_dofs(),
600 ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
601
602 std::vector<VectorType> all_in = {in};
603 std::vector<VectorType> all_out = {out};
604
605 interpolate(all_in, all_out);
606
607 out = all_out[0];
608}
609
610
611
612template <int dim, typename VectorType, int spacedim>
613std::size_t
615{
616 // at the moment we do not include the memory
617 // consumption of the cell_map as we have no
618 // real idea about memory consumption of a
619 // std::map
620 return (MemoryConsumption::memory_consumption(dof_handler) +
622 sizeof(prepared_for) +
624 MemoryConsumption::memory_consumption(dof_values_on_cell));
625}
626
627
628
629template <int dim, typename VectorType, int spacedim>
630std::size_t
636
637
638/*-------------- Explicit Instantiations -------------------------------*/
639#define SPLIT_INSTANTIATIONS_COUNT 4
640#ifndef SPLIT_INSTANTIATIONS_INDEX
641# define SPLIT_INSTANTIATIONS_INDEX 0
642#endif
643#include "solution_transfer.inst"
644
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
SmartPointer< const DoFHandler< dim, spacedim >, SolutionTransfer< dim, VectorType, spacedim > > dof_handler
std::size_t memory_consumption() const
void prepare_for_coarsening_and_refinement(const std::vector< VectorType > &all_in)
SolutionTransfer(const DoFHandler< dim, spacedim > &dof)
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
void swap(SmartPointer< T, P > &t1, SmartPointer< T, Q > &t2)
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)