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