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