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