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