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