deal.II version GIT relicensing-2850-g1646a780ac 2025-03-17 15: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
dof_accessor_get.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2023 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
18
19#include <deal.II/fe/fe.h>
20
22
33#include <deal.II/lac/vector.h>
34
35#include <vector>
36
38
39
40template <int dim, int spacedim, bool lda>
41template <typename Number>
42void
44 const ReadVector<Number> &values,
45 Vector<Number> &interpolated_values,
46 const types::fe_index fe_index_) const
47{
48 const types::fe_index fe_index =
49 (this->dof_handler->hp_capability_enabled == false &&
50 fe_index_ == numbers::invalid_fe_index) ?
52 fe_index_;
53
54 if (this->is_active())
55 // If this cell is active: simply return the exact values on this
56 // cell unless the finite element we need to interpolate to is different
57 // than the one we have on the current cell
58 {
59 if ((this->dof_handler->hp_capability_enabled == false) ||
60 // for hp-DoFHandlers, we need to require that on
61 // active cells, you either don't specify an fe_index,
62 // or that you specify the correct one
63 (fe_index == this->active_fe_index()) ||
64 (fe_index == numbers::invalid_fe_index))
65 this->get_dof_values(values, interpolated_values);
66 else
67 {
68 // well, here we need to first get the values from the current
69 // cell and then interpolate it to the element requested. this
70 // can clearly only happen for DoFHandler objects in hp-mode
71 const unsigned int dofs_per_cell = this->get_fe().n_dofs_per_cell();
72 if (dofs_per_cell == 0)
73 {
74 interpolated_values = 0;
75 }
76 else
77 {
78 Vector<Number> tmp(dofs_per_cell);
79 this->get_dof_values(values, tmp);
80
81 FullMatrix<double> interpolation(
82 this->dof_handler->get_fe(fe_index).n_dofs_per_cell(),
83 this->get_fe().n_dofs_per_cell());
84 this->dof_handler->get_fe(fe_index).get_interpolation_matrix(
85 this->get_fe(), interpolation);
86 interpolation.vmult(interpolated_values, tmp);
87 }
88 }
89 }
90 else
91 // The cell is not active; we need to obtain data them from
92 // children recursively.
93 {
94 // we are on a non-active cell. these do not have any finite
95 // element associated with them in the hp-context (in the non-hp-
96 // context, we can simply assume that the FE space to which we
97 // want to interpolate is the same as for all elements in the
98 // mesh). consequently, we cannot interpolate from children's FE
99 // space to this cell's (unknown) FE space unless an explicit
100 // fe_index is given
101 Assert((this->dof_handler->hp_capability_enabled == false) ||
102 (fe_index != numbers::invalid_fe_index),
104 "You cannot call this function on non-active cells "
105 "of DoFHandler objects unless you provide an explicit "
106 "finite element index because they do not have naturally "
107 "associated finite element spaces associated: degrees "
108 "of freedom are only distributed on active cells for which "
109 "the active FE index has been set."));
110
112 this->get_dof_handler().get_fe(fe_index);
113 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
114
115 Assert(this->dof_handler != nullptr,
116 typename BaseClass::ExcInvalidObject());
117 Assert(interpolated_values.size() == dofs_per_cell,
118 typename BaseClass::ExcVectorDoesNotMatch());
119 Assert(values.size() == this->dof_handler->n_dofs(),
120 typename BaseClass::ExcVectorDoesNotMatch());
121
122
123 // see if the finite element we have on the current cell has any
124 // degrees of freedom to begin with; if not (e.g., when
125 // interpolating FE_Nothing), then simply skip all of the
126 // following since the output vector would be of size zero
127 // anyway (and in fact is of size zero, see the assertion above)
128 if (fe.n_dofs_per_cell() > 0)
129 {
130 Vector<Number> tmp1(dofs_per_cell);
131 Vector<Number> tmp2(dofs_per_cell);
132
133 interpolated_values = 0;
134
135 // later on we will have to push the values interpolated from the
136 // child to the mother cell into the output vector. unfortunately,
137 // there are two types of elements: ones where you add up the
138 // contributions from the different child cells, and ones where you
139 // overwrite.
140 //
141 // an example for the first is piecewise constant (and discontinuous)
142 // elements, where we build the value on the coarse cell by averaging
143 // the values from the cell (i.e. by adding up a fraction of the
144 // values of their values)
145 //
146 // an example for the latter are the usual continuous elements. the
147 // value on a vertex of a coarse cell must there be the same,
148 // irrespective of the adjacent cell we are presently on. so we always
149 // overwrite. in fact, we must, since we cannot know in advance how
150 // many neighbors there will be, so there is no way to compute the
151 // average with fixed factors
152 //
153 // so we have to find out to which type this element belongs. the
154 // difficulty is: the finite element may be a composed one, so we can
155 // only hope to do this for each shape function individually. in fact,
156 // there are even weird finite elements (for example the
157 // Raviart-Thomas element) which have shape functions that are
158 // additive (interior ones) and others that are overwriting (face
159 // degrees of freedom that need to be continuous across the face).
160 for (unsigned int child = 0; child < this->n_children(); ++child)
161 {
162 // get the values from the present child, if necessary by
163 // interpolation itself either from its own children or
164 // by interpolating from the finite element on an active
165 // child to the finite element space requested here
166 this->child(child)->get_interpolated_dof_values(values,
167 tmp1,
168 fe_index);
169 // interpolate these to the mother cell
170 fe.get_restriction_matrix(child, this->refinement_case())
171 .vmult(tmp2, tmp1);
172
173 // and add up or set them in the output vector
174 for (unsigned int i = 0; i < dofs_per_cell; ++i)
175 if (fe.restriction_is_additive(i))
176 interpolated_values(i) += tmp2(i);
177 else if (tmp2(i) != Number())
178 interpolated_values(i) = tmp2(i);
179 }
180 }
181 }
182}
183
184
185// --------------------------------------------------------------------------
186// explicit instantiations
187#include "dofs/dof_accessor_get.inst"
188
void get_interpolated_dof_values(const ReadVector< Number > &values, Vector< Number > &interpolated_values, const types::fe_index fe_index=numbers::invalid_fe_index) const
unsigned int n_dofs_per_cell() const
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
bool restriction_is_additive(const unsigned int index) const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr types::fe_index invalid_fe_index
Definition types.h:254
unsigned short int fe_index
Definition types.h:68