Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18: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 <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<Number>::value,
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<Number>::value, 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<VectorType,
81  std::is_same<VectorType,
83  std::is_same<
84  VectorType,
86  std::is_same<VectorType,
88  typename VectorType::value_type>>::value ||
89  std::is_same<VectorType,
91  typename VectorType::value_type>>::value;
92 
97  template <typename T>
99  decltype(std::declval<T const>().set_ghost_state(std::declval<bool>()));
100 
101  template <typename T>
102  constexpr bool has_set_ghost_state =
103  is_supported_operation<set_ghost_state_t, T>;
104 
105  template <typename VectorType,
106  typename std::enable_if<has_set_ghost_state<VectorType>,
107  VectorType>::type * = nullptr>
108  void
109  set_ghost_state(VectorType &vector, const bool ghosted)
110  {
111  vector.set_ghost_state(ghosted);
112  }
113 
114  template <typename VectorType,
115  typename std::enable_if<!has_set_ghost_state<VectorType>,
116  VectorType>::type * = nullptr>
117  void
118  set_ghost_state(VectorType &, const bool)
119  {
120  // serial vector: nothing to do
121  }
122 #endif
123 
128  template <int dim,
129  int spacedim,
130  bool lda,
131  class OutputVector,
132  typename number>
133  void
135  const Vector<number> & local_values,
136  OutputVector & values,
137  const bool perform_check)
138  {
139  (void)perform_check;
140 
141 #ifdef DEBUG
142  if (perform_check && is_dealii_vector<OutputVector>)
143  {
144  const bool old_ghost_state = values.has_ghost_elements();
145  set_ghost_state(values, true);
146 
147  Vector<number> local_values_old(cell.get_fe().n_dofs_per_cell());
148  cell.get_dof_values(values, local_values_old);
149 
150  for (unsigned int i = 0; i < cell.get_fe().n_dofs_per_cell(); ++i)
151  {
152  // a check consistent with the one in
153  // Utilities::MPI::Partitioner::import_from_ghosted_array_finish()
154  Assert(local_values_old[i] == number() ||
155  get_abs(local_values_old[i] - local_values[i]) <=
156  get_abs(local_values_old[i] + local_values[i]) *
157  100000. *
158  std::numeric_limits<typename numbers::NumberTraits<
159  number>::real_type>::epsilon(),
160  ExcNonMatchingElementsSetDofValuesByInterpolation<number>(
161  local_values[i], local_values_old[i]));
162  }
163 
164  set_ghost_state(values, old_ghost_state);
165  }
166 #endif
167 
168  cell.set_dof_values(local_values, values);
169  }
170 
171 
172  template <int dim,
173  int spacedim,
174  bool lda,
175  class OutputVector,
176  typename number>
177  void
180  const Vector<number> & local_values,
181  OutputVector & values,
182  const unsigned int fe_index_,
183  const std::function<void(const DoFCellAccessor<dim, spacedim, lda> &cell,
184  const Vector<number> &local_values,
185  OutputVector & values)> &processor)
186  {
187  const unsigned int fe_index =
188  (cell.get_dof_handler().has_hp_capabilities() == false &&
191  fe_index_;
192 
193  if (cell.is_active() && !cell.is_artificial())
194  {
195  if ((cell.get_dof_handler().has_hp_capabilities() == false) ||
196  // for hp-DoFHandlers, we need to require that on
197  // active cells, you either don't specify an fe_index,
198  // or that you specify the correct one
199  (fe_index == cell.active_fe_index()) ||
201  // simply set the values on this cell
202  processor(cell, local_values, values);
203  else
204  {
205  Assert(local_values.size() ==
206  cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell(),
207  ExcMessage("Incorrect size of local_values vector."));
208 
209  FullMatrix<double> interpolation(
210  cell.get_fe().n_dofs_per_cell(),
211  cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell());
212 
213  cell.get_fe().get_interpolation_matrix(
214  cell.get_dof_handler().get_fe(fe_index), interpolation);
215 
216  // do the interpolation to the target space. for historical
217  // reasons, matrices are set to size 0x0 internally even if
218  // we reinit as 4x0, so we have to treat this case specially
219  Vector<number> tmp(cell.get_fe().n_dofs_per_cell());
220  if ((tmp.size() > 0) && (local_values.size() > 0))
221  interpolation.vmult(tmp, local_values);
222 
223  // now set the dof values in the global vector
224  processor(cell, tmp, values);
225  }
226  }
227  else
228  // otherwise distribute them to the children
229  {
230  Assert((cell.get_dof_handler().has_hp_capabilities() == false) ||
232  ExcMessage(
233  "You cannot call this function on non-active cells "
234  "of DoFHandler objects unless you provide an explicit "
235  "finite element index because they do not have naturally "
236  "associated finite element spaces associated: degrees "
237  "of freedom are only distributed on active cells for which "
238  "the active FE index has been set."));
239 
240  const FiniteElement<dim, spacedim> &fe =
241  cell.get_dof_handler().get_fe(fe_index);
242  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
243 
244  Assert(local_values.size() == dofs_per_cell,
246  ExcVectorDoesNotMatch()));
247  Assert(values.size() == cell.get_dof_handler().n_dofs(),
249  ExcVectorDoesNotMatch()));
250 
251  Vector<number> tmp(dofs_per_cell);
252 
253  for (unsigned int child = 0; child < cell.n_children(); ++child)
254  {
255  if (tmp.size() > 0)
256  fe.get_prolongation_matrix(child, cell.refinement_case())
257  .vmult(tmp, local_values);
259  *cell.child(child), tmp, values, fe_index, processor);
260  }
261  }
262  }
263 
264 } // namespace internal
265 
266 
267 
268 template <int dim, int spacedim, bool lda>
269 template <class OutputVector, typename number>
270 void
272  const Vector<number> &local_values,
273  OutputVector & values,
274  const unsigned int fe_index_,
275  const bool perform_check) const
276 {
277  internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
278  *this,
279  local_values,
280  values,
281  fe_index_,
282  [perform_check](const DoFCellAccessor<dim, spacedim, lda> &cell,
283  const Vector<number> & local_values,
284  OutputVector & values) {
285  internal::set_dof_values(cell, local_values, values, perform_check);
286  });
287 }
288 
289 
290 template <int dim, int spacedim, bool lda>
291 template <class OutputVector, typename number>
292 void
295  const Vector<number> &local_values,
296  OutputVector & values,
297  const unsigned int fe_index_) const
298 {
299  internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
300  *this,
301  local_values,
302  values,
303  fe_index_,
305  const Vector<number> & local_values,
306  OutputVector & values) {
307  std::vector<types::global_dof_index> dof_indices(
308  cell.get_fe().n_dofs_per_cell());
309  cell.get_dof_indices(dof_indices);
311  dof_indices,
312  values);
313  });
314 }
315 
316 
317 // --------------------------------------------------------------------------
318 // explicit instantiations
319 #include "dof_accessor_set.inst"
320 
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
void distribute_local_to_global_by_interpolation(const Vector< number > &local_values, OutputVector &values, const unsigned int fe_index=DoFHandler< dimension_, space_dimension_ >::invalid_fe_index) const
unsigned int active_fe_index() const
TriaIterator< DoFCellAccessor< dimension_, space_dimension_, level_dof_access > > child(const unsigned int i) const
void set_dof_values(const Vector< number > &local_values, OutputVector &values) const
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
void set_dof_values_by_interpolation(const Vector< number > &local_values, OutputVector &values, const unsigned int fe_index=DoFHandler< dimension_, space_dimension_ >::invalid_fe_index, const bool perform_check=false) const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int 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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNonMatchingElementsSetDofValuesByInterpolation(Number arg1, Number arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type size() const
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
void process_by_interpolation(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values, const unsigned int fe_index_, const std::function< void(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values)> &processor)
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)