Reference documentation for deal.II version GIT 20f059c89a 2022-12-04 14:55: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_set.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 
28 #include <deal.II/lac/la_vector.h>
36 #include <deal.II/lac/vector.h>
37 
38 #include <limits>
39 #include <vector>
40 
42 
43 template <typename Number>
45  Number,
46  Number,
47  << "Called set_dof_values_by_interpolation(), but"
48  << " the element to be set, value " << std::setprecision(16)
49  << arg1 << ", does not match with the non-zero value "
50  << std::setprecision(16) << arg2 << " already set before.");
51 
52 namespace internal
53 {
54 #ifdef DEBUG
60  template <typename Number>
61  std::enable_if_t<!std::is_unsigned<Number>::value,
63  get_abs(const Number a)
64  {
65  return std::abs(a);
66  }
67 
68  template <typename Number>
69  std::enable_if_t<std::is_unsigned<Number>::value, Number>
70  get_abs(const Number a)
71  {
72  return a;
73  }
74 
78  template <typename VectorType>
79  constexpr bool is_dealii_vector =
80  std::is_same<VectorType,
82  std::is_same<VectorType,
84  std::is_same<
85  VectorType,
87  std::is_same<VectorType,
89  typename VectorType::value_type>>::value ||
90  std::is_same<VectorType,
92  typename VectorType::value_type>>::value;
93 
98  template <typename T>
100  decltype(std::declval<T const>().set_ghost_state(std::declval<bool>()));
101 
102  template <typename T>
103  constexpr bool has_set_ghost_state =
104  is_supported_operation<set_ghost_state_t, T>;
105 
106  template <
107  typename VectorType,
108  std::enable_if_t<has_set_ghost_state<VectorType>, VectorType> * = nullptr>
109  void
110  set_ghost_state(VectorType &vector, const bool ghosted)
111  {
112  vector.set_ghost_state(ghosted);
113  }
114 
115  template <
116  typename VectorType,
117  std::enable_if_t<!has_set_ghost_state<VectorType>, VectorType> * = nullptr>
118  void
119  set_ghost_state(VectorType &, const bool)
120  {
121  // serial vector: nothing to do
122  }
123 #endif
124 
129  template <int dim,
130  int spacedim,
131  bool lda,
132  class OutputVector,
133  typename number>
134  void
136  const Vector<number> & local_values,
137  OutputVector & values,
138  const bool perform_check)
139  {
140  (void)perform_check;
141 
142 #ifdef DEBUG
143  if (perform_check && is_dealii_vector<OutputVector>)
144  {
145  const bool old_ghost_state = values.has_ghost_elements();
146  set_ghost_state(values, true);
147 
148  Vector<number> local_values_old(cell.get_fe().n_dofs_per_cell());
149  cell.get_dof_values(values, local_values_old);
150 
151  for (unsigned int i = 0; i < cell.get_fe().n_dofs_per_cell(); ++i)
152  {
153  // a check consistent with the one in
154  // Utilities::MPI::Partitioner::import_from_ghosted_array_finish()
155  Assert(local_values_old[i] == number() ||
156  get_abs(local_values_old[i] - local_values[i]) <=
157  get_abs(local_values_old[i] + local_values[i]) *
158  100000. *
159  std::numeric_limits<typename numbers::NumberTraits<
160  number>::real_type>::epsilon(),
161  ExcNonMatchingElementsSetDofValuesByInterpolation<number>(
162  local_values[i], local_values_old[i]));
163  }
164 
165  set_ghost_state(values, old_ghost_state);
166  }
167 #endif
168 
169  cell.set_dof_values(local_values, values);
170  }
171 
172 
173  template <int dim,
174  int spacedim,
175  bool lda,
176  class OutputVector,
177  typename number>
178  void
181  const Vector<number> & local_values,
182  OutputVector & values,
183  const types::fe_index fe_index_,
184  const std::function<void(const DoFCellAccessor<dim, spacedim, lda> &cell,
185  const Vector<number> &local_values,
186  OutputVector & values)> &processor)
187  {
188  const types::fe_index fe_index =
189  (cell.get_dof_handler().has_hp_capabilities() == false &&
190  fe_index_ == numbers::invalid_fe_index) ?
192  fe_index_;
193 
194  if (cell.is_active() && !cell.is_artificial())
195  {
196  if ((cell.get_dof_handler().has_hp_capabilities() == false) ||
197  // for hp-DoFHandlers, we need to require that on
198  // active cells, you either don't specify an fe_index,
199  // or that you specify the correct one
200  (fe_index == cell.active_fe_index()) ||
202  // simply set the values on this cell
203  processor(cell, local_values, values);
204  else
205  {
206  Assert(local_values.size() ==
208  ExcMessage("Incorrect size of local_values vector."));
209 
210  FullMatrix<double> interpolation(
211  cell.get_fe().n_dofs_per_cell(),
213 
214  cell.get_fe().get_interpolation_matrix(
215  cell.get_dof_handler().get_fe(fe_index), interpolation);
216 
217  // do the interpolation to the target space. for historical
218  // reasons, matrices are set to size 0x0 internally even if
219  // we reinit as 4x0, so we have to treat this case specially
220  Vector<number> tmp(cell.get_fe().n_dofs_per_cell());
221  if ((tmp.size() > 0) && (local_values.size() > 0))
222  interpolation.vmult(tmp, local_values);
223 
224  // now set the dof values in the global vector
225  processor(cell, tmp, values);
226  }
227  }
228  else
229  // otherwise distribute them to the children
230  {
231  Assert((cell.get_dof_handler().has_hp_capabilities() == false) ||
233  ExcMessage(
234  "You cannot call this function on non-active cells "
235  "of DoFHandler objects unless you provide an explicit "
236  "finite element index because they do not have naturally "
237  "associated finite element spaces associated: degrees "
238  "of freedom are only distributed on active cells for which "
239  "the active FE index has been set."));
240 
241  const FiniteElement<dim, spacedim> &fe =
243  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
244 
245  Assert(local_values.size() == dofs_per_cell,
247  ExcVectorDoesNotMatch()));
248  Assert(values.size() == cell.get_dof_handler().n_dofs(),
250  ExcVectorDoesNotMatch()));
251 
252  Vector<number> tmp(dofs_per_cell);
253 
254  for (unsigned int child = 0; child < cell.n_children(); ++child)
255  {
256  if (tmp.size() > 0)
257  fe.get_prolongation_matrix(child, cell.refinement_case())
258  .vmult(tmp, local_values);
260  *cell.child(child), tmp, values, fe_index, processor);
261  }
262  }
263  }
264 
265 } // namespace internal
266 
267 
268 
269 template <int dim, int spacedim, bool lda>
270 template <class OutputVector, typename number>
271 void
273  const Vector<number> &local_values,
274  OutputVector & values,
275  const types::fe_index fe_index_,
276  const bool perform_check) const
277 {
278  internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
279  *this,
280  local_values,
281  values,
282  fe_index_,
283  [perform_check](const DoFCellAccessor<dim, spacedim, lda> &cell,
284  const Vector<number> & local_values,
285  OutputVector & values) {
286  internal::set_dof_values(cell, local_values, values, perform_check);
287  });
288 }
289 
290 
291 template <int dim, int spacedim, bool lda>
292 template <class OutputVector, typename number>
293 void
296  const Vector<number> &local_values,
297  OutputVector & values,
298  const types::fe_index fe_index_) const
299 {
300  internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
301  *this,
302  local_values,
303  values,
304  fe_index_,
306  const Vector<number> & local_values,
307  OutputVector & values) {
308  std::vector<types::global_dof_index> dof_indices(
309  cell.get_fe().n_dofs_per_cell());
310  cell.get_dof_indices(dof_indices);
312  dof_indices,
313  values);
314  });
315 }
316 
317 
318 // --------------------------------------------------------------------------
319 // explicit instantiations
320 #include "dof_accessor_set.inst"
321 
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
const DoFHandler< dim, spacedim > & get_dof_handler() const
void get_dof_values(const InputVector &values, Vector< number > &local_values) const
TriaIterator< DoFCellAccessor< dimension_, space_dimension_, level_dof_access > > child(const unsigned int i) const
void distribute_local_to_global_by_interpolation(const Vector< number > &local_values, OutputVector &values, const types::fe_index fe_index=numbers::invalid_fe_index) const
void set_dof_values(const Vector< number > &local_values, OutputVector &values) const
void set_dof_values_by_interpolation(const Vector< number > &local_values, OutputVector &values, const types::fe_index fe_index=numbers::invalid_fe_index, const bool perform_check=false) const
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
types::fe_index active_fe_index() const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
unsigned int n_dofs_per_cell() const
virtual const FullMatrix< double > & get_prolongation_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
size_type size() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcNonMatchingElementsSetDofValuesByInterpolation(Number arg1, Number arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
std::enable_if_t<!std::is_unsigned< Number >::value, typename numbers::NumberTraits< Number >::real_type > get_abs(const Number a)
constexpr bool has_set_ghost_state
decltype(std::declval< T const >().set_ghost_state(std::declval< bool >())) set_ghost_state_t
constexpr bool is_dealii_vector
void set_ghost_state(VectorType &vector, const bool ghosted)
void set_dof_values(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values, const bool perform_check)
void process_by_interpolation(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values, const types::fe_index fe_index_, const std::function< void(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values)> &processor)
const types::fe_index invalid_fe_index
Definition: types.h:221
unsigned short int fe_index
Definition: types.h:59