Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
36#include <deal.II/lac/vector.h>
37
38#include <limits>
39#include <vector>
40
42
43template <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
52namespace 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()) ||
201 (fe_index == numbers::invalid_fe_index))
202 // simply set the values on this cell
203 processor(cell, local_values, values);
204 else
205 {
206 Assert(local_values.size() ==
207 cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell(),
208 ExcMessage("Incorrect size of local_values vector."));
209
210 FullMatrix<double> interpolation(
211 cell.get_fe().n_dofs_per_cell(),
212 cell.get_dof_handler().get_fe(fe_index).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) ||
232 (fe_index != numbers::invalid_fe_index),
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
242 cell.get_dof_handler().get_fe(fe_index);
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
269template <int dim, int spacedim, bool lda>
270template <class OutputVector, typename number>
271void
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
291template <int dim, int spacedim, bool lda>
292template <class OutputVector, typename number>
293void
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
const FiniteElement< dimension_, space_dimension_ > & get_fe() 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
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
types::fe_index active_fe_index() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
static ::ExceptionBase & ExcNonMatchingElementsSetDofValuesByInterpolation(Number arg1, Number arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
std::enable_if_t<!std::is_unsigned< Number >::value, typename numbers::NumberTraits< Number >::real_type > get_abs(const Number a)
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:228
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:60