Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
dof_accessor_get.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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 
16 #include <deal.II/dofs/dof_accessor.h>
17 #include <deal.II/dofs/dof_handler.h>
18 #include <deal.II/dofs/dof_levels.h>
19 
20 #include <deal.II/fe/fe.h>
21 
22 #include <deal.II/grid/tria_iterator.h>
23 #include <deal.II/grid/tria_iterator.templates.h>
24 
25 #include <deal.II/hp/dof_handler.h>
26 
27 #include <deal.II/lac/block_vector.h>
28 #include <deal.II/lac/la_parallel_block_vector.h>
29 #include <deal.II/lac/la_parallel_vector.h>
30 #include <deal.II/lac/la_vector.h>
31 #include <deal.II/lac/petsc_block_vector.h>
32 #include <deal.II/lac/petsc_vector.h>
33 #include <deal.II/lac/sparse_matrix.h>
34 #include <deal.II/lac/trilinos_parallel_block_vector.h>
35 #include <deal.II/lac/trilinos_tpetra_vector.h>
36 #include <deal.II/lac/trilinos_vector.h>
37 #include <deal.II/lac/vector.h>
38 
39 #include <vector>
40 
41 DEAL_II_NAMESPACE_OPEN
42 
43 
44 template <typename DoFHandlerType, bool lda>
45 template <class InputVector, typename number>
46 void
48  const InputVector &values,
49  Vector<number> & interpolated_values,
50  const unsigned int fe_index) const
51 {
52  if (!this->has_children())
53  // if this cell has no children: simply return the exact values on this
54  // cell unless the finite element we need to interpolate to is different
55  // than the one we have on the current cell
56  {
57  if ((dynamic_cast<DoFHandler<DoFHandlerType::dimension,
58  DoFHandlerType::space_dimension> *>(
59  this->dof_handler) != nullptr) ||
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 == DoFHandlerType::default_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 hp::DoFHandler objects
71  Vector<number> tmp(this->get_fe().dofs_per_cell);
72  this->get_dof_values(values, tmp);
73 
74  FullMatrix<double> interpolation(
75  this->dof_handler->get_fe(fe_index).dofs_per_cell,
76  this->get_fe().dofs_per_cell);
77  this->dof_handler->get_fe(fe_index).get_interpolation_matrix(
78  this->get_fe(), interpolation);
79  interpolation.vmult(interpolated_values, tmp);
80  }
81  }
82  else
83  // otherwise obtain them from the children
84  {
85  // we are on a non-active cell. these do not have any finite
86  // element associated with them in the hp context (in the non-hp
87  // context, we can simply assume that the FE space to which we
88  // want to interpolate is the same as for all elements in the
89  // mesh). consequently, we cannot interpolate from children's FE
90  // space to this cell's (unknown) FE space unless an explicit
91  // fe_index is given
92  Assert((dynamic_cast<DoFHandler<DoFHandlerType::dimension,
93  DoFHandlerType::space_dimension> *>(
94  this->dof_handler) != nullptr) ||
95  (fe_index != DoFHandlerType::default_fe_index),
96  ExcMessage(
97  "You cannot call this function on non-active cells "
98  "of hp::DoFHandler objects unless you provide an explicit "
99  "finite element index because they do not have naturally "
100  "associated finite element spaces associated: degrees "
101  "of freedom are only distributed on active cells for which "
102  "the active_fe_index has been set."));
103 
104  const FiniteElement<dim, spacedim> &fe =
105  this->get_dof_handler().get_fe(fe_index);
106  const unsigned int dofs_per_cell = fe.dofs_per_cell;
107 
108  Assert(this->dof_handler != nullptr,
109  typename BaseClass::ExcInvalidObject());
110  Assert(interpolated_values.size() == dofs_per_cell,
111  typename BaseClass::ExcVectorDoesNotMatch());
112  Assert(values.size() == this->dof_handler->n_dofs(),
113  typename BaseClass::ExcVectorDoesNotMatch());
114 
115 
116  // see if the finite element we have on the current cell has any
117  // degrees of freedom to begin with; if not (e.g., when
118  // interpolating FE_Nothing), then simply skip all of the
119  // following since the output vector would be of size zero
120  // anyway (and in fact is of size zero, see the assertion above)
121  if (fe.dofs_per_cell > 0)
122  {
123  Vector<number> tmp1(dofs_per_cell);
124  Vector<number> tmp2(dofs_per_cell);
125 
126  interpolated_values = 0;
127 
128  // later on we will have to push the values interpolated from the
129  // child to the mother cell into the output vector. unfortunately,
130  // there are two types of elements: ones where you add up the
131  // contributions from the different child cells, and ones where you
132  // overwrite.
133  //
134  // an example for the first is piecewise constant (and discontinuous)
135  // elements, where we build the value on the coarse cell by averaging
136  // the values from the cell (i.e. by adding up a fraction of the
137  // values of their values)
138  //
139  // an example for the latter are the usual continuous elements. the
140  // value on a vertex of a coarse cell must there be the same,
141  // irrespective of the adjacent cell we are presently on. so we always
142  // overwrite. in fact, we must, since we cannot know in advance how
143  // many neighbors there will be, so there is no way to compute the
144  // average with fixed factors
145  //
146  // so we have to find out to which type this element belongs. the
147  // difficulty is: the finite element may be a composed one, so we can
148  // only hope to do this for each shape function individually. in fact,
149  // there are even weird finite elements (for example the
150  // Raviart-Thomas element) which have shape functions that are
151  // additive (interior ones) and others that are overwriting (face
152  // degrees of freedom that need to be continuous across the face).
153  for (unsigned int child = 0; child < this->n_children(); ++child)
154  {
155  // get the values from the present child, if necessary by
156  // interpolation itself either from its own children or
157  // by interpolating from the finite element on an active
158  // child to the finite element space requested here
159  this->child(child)->get_interpolated_dof_values(values,
160  tmp1,
161  fe_index);
162  // interpolate these to the mother cell
163  fe.get_restriction_matrix(child, this->refinement_case())
164  .vmult(tmp2, tmp1);
165 
166  // and add up or set them in the output vector
167  for (unsigned int i = 0; i < dofs_per_cell; ++i)
168  if (fe.restriction_is_additive(i))
169  interpolated_values(i) += tmp2(i);
170  else if (tmp2(i) != number())
171  interpolated_values(i) = tmp2(i);
172  }
173  }
174  }
175 }
176 
177 
178 // --------------------------------------------------------------------------
179 // explicit instantiations
180 #include "dof_accessor_get.inst"
181 
182 DEAL_II_NAMESPACE_CLOSE
bool restriction_is_additive(const unsigned int index) const
Definition: fe.h:3276
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void get_interpolated_dof_values(const InputVector &values, Vector< number > &interpolated_values, const unsigned int fe_index=DoFHandlerType::default_fe_index) const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:306
const unsigned int dofs_per_cell
Definition: fe_base.h:282
size_type size() const