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