deal.II version GIT relicensing-2850-g1646a780ac 2025-03-17 15:10:00+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\}}\)
Loading...
Searching...
No Matches
dof_accessor_set.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2013 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
18
19#include <deal.II/fe/fe.h>
20
22
33#include <deal.II/lac/vector.h>
34
35#include <iomanip>
36#include <limits>
37#include <vector>
38
40
41template <typename Number>
43 Number,
44 Number,
45 << "Called set_dof_values_by_interpolation(), but"
46 << " the element to be set, value " << std::setprecision(16)
47 << arg1 << ", does not match with the non-zero value "
48 << std::setprecision(16) << arg2 << " already set before.");
49
50namespace internal
51{
52#ifdef DEBUG
58 template <typename Number>
59 std::enable_if_t<!std::is_unsigned_v<Number>,
61 get_abs(const Number a)
62 {
63 return std::abs(a);
64 }
65
66 template <typename Number>
67 std::enable_if_t<std::is_unsigned_v<Number>, Number>
68 get_abs(const Number a)
69 {
70 return a;
71 }
72
73
78 template <typename T>
80 decltype(std::declval<const T>().set_ghost_state(std::declval<bool>()));
81
82 template <typename T>
83 constexpr bool has_set_ghost_state =
84 is_supported_operation<set_ghost_state_t, T>;
85
86 template <
87 typename VectorType,
88 std::enable_if_t<has_set_ghost_state<VectorType>, VectorType> * = nullptr>
89 void
90 set_ghost_state(VectorType &vector, const bool ghosted)
91 {
92 vector.set_ghost_state(ghosted);
93 }
94
95 template <
96 typename VectorType,
97 std::enable_if_t<!has_set_ghost_state<VectorType>, VectorType> * = nullptr>
98 void
99 set_ghost_state(VectorType &, const bool)
100 {
101 // serial vector: nothing to do
102 }
103#endif
104
105 namespace
106 {
107 // Test whether a vector is a deal.II vector
108 template <typename VectorType>
109 constexpr bool is_dealii_vector =
110 std::is_same_v<VectorType,
112 std::is_same_v<VectorType,
114 std::is_same_v<VectorType,
116 typename VectorType::value_type>> ||
117 std::is_same_v<VectorType,
119 typename VectorType::value_type>>;
120 } // namespace
121
122
127 template <int dim,
128 int spacedim,
129 bool lda,
130 class OutputVector,
131 typename number>
132 void
134 const Vector<number> &local_values,
135 OutputVector &values,
136 const bool perform_check)
137 {
138 (void)perform_check;
139
140 if constexpr (running_in_debug_mode())
141 {
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(
155 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]) * 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 }
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 types::fe_index 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 types::fe_index fe_index =
188 (cell.get_dof_handler().has_hp_capabilities() == false &&
189 fe_index_ == numbers::invalid_fe_index) ?
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()) ||
200 (fe_index == numbers::invalid_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) ||
231 (fe_index != numbers::invalid_fe_index),
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
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
268template <int dim, int spacedim, bool lda>
269template <class OutputVector, typename number>
270void
272 const Vector<number> &local_values,
273 OutputVector &values,
274 const types::fe_index 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
290template <int dim, int spacedim, bool lda>
291template <class OutputVector, typename number>
292void
295 const Vector<number> &local_values,
296 OutputVector &values,
297 const types::fe_index 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 "dofs/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
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
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ExcNonMatchingElementsSetDofValuesByInterpolation(Number arg1, Number arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
decltype(std::declval< const T >().set_ghost_state(std::declval< bool >())) set_ghost_state_t
constexpr bool has_set_ghost_state
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)
constexpr types::fe_index invalid_fe_index
Definition types.h:254
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:68