deal.II version GIT relicensing-1927-g3de9220933 2024-10-03 08:40: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#include <deal.II/grid/tria_iterator.templates.h>
23
34#include <deal.II/lac/vector.h>
35
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
76 template <typename VectorType>
77 constexpr bool is_dealii_vector =
78 std::is_same_v<VectorType,
80 std::is_same_v<VectorType,
82 std::is_same_v<VectorType,
84 typename VectorType::value_type>> ||
85 std::is_same_v<VectorType,
87 typename VectorType::value_type>>;
88
93 template <typename T>
95 decltype(std::declval<const T>().set_ghost_state(std::declval<bool>()));
96
97 template <typename T>
98 constexpr bool has_set_ghost_state =
99 is_supported_operation<set_ghost_state_t, T>;
100
101 template <
102 typename VectorType,
103 std::enable_if_t<has_set_ghost_state<VectorType>, VectorType> * = nullptr>
104 void
105 set_ghost_state(VectorType &vector, const bool ghosted)
106 {
107 vector.set_ghost_state(ghosted);
108 }
109
110 template <
111 typename VectorType,
112 std::enable_if_t<!has_set_ghost_state<VectorType>, VectorType> * = nullptr>
113 void
114 set_ghost_state(VectorType &, const bool)
115 {
116 // serial vector: nothing to do
117 }
118#endif
119
124 template <int dim,
125 int spacedim,
126 bool lda,
127 class OutputVector,
128 typename number>
129 void
131 const Vector<number> &local_values,
132 OutputVector &values,
133 const bool perform_check)
134 {
135 (void)perform_check;
136
137#ifdef DEBUG
138 if (perform_check && is_dealii_vector<OutputVector>)
139 {
140 const bool old_ghost_state = values.has_ghost_elements();
141 set_ghost_state(values, true);
142
143 Vector<number> local_values_old(cell.get_fe().n_dofs_per_cell());
144 cell.get_dof_values(values, local_values_old);
145
146 for (unsigned int i = 0; i < cell.get_fe().n_dofs_per_cell(); ++i)
147 {
148 // a check consistent with the one in
149 // Utilities::MPI::Partitioner::import_from_ghosted_array_finish()
150 Assert(local_values_old[i] == number() ||
151 get_abs(local_values_old[i] - local_values[i]) <=
152 get_abs(local_values_old[i] + local_values[i]) *
153 100000. *
154 std::numeric_limits<typename numbers::NumberTraits<
155 number>::real_type>::epsilon(),
156 ExcNonMatchingElementsSetDofValuesByInterpolation<number>(
157 local_values[i], local_values_old[i]));
158 }
159
160 set_ghost_state(values, old_ghost_state);
161 }
162#endif
163
164 cell.set_dof_values(local_values, values);
165 }
166
167
168 template <int dim,
169 int spacedim,
170 bool lda,
171 class OutputVector,
172 typename number>
173 void
176 const Vector<number> &local_values,
177 OutputVector &values,
178 const types::fe_index fe_index_,
179 const std::function<void(const DoFCellAccessor<dim, spacedim, lda> &cell,
180 const Vector<number> &local_values,
181 OutputVector &values)> &processor)
182 {
183 const types::fe_index fe_index =
184 (cell.get_dof_handler().has_hp_capabilities() == false &&
185 fe_index_ == numbers::invalid_fe_index) ?
187 fe_index_;
188
189 if (cell.is_active() && !cell.is_artificial())
190 {
191 if ((cell.get_dof_handler().has_hp_capabilities() == false) ||
192 // for hp-DoFHandlers, we need to require that on
193 // active cells, you either don't specify an fe_index,
194 // or that you specify the correct one
195 (fe_index == cell.active_fe_index()) ||
196 (fe_index == numbers::invalid_fe_index))
197 // simply set the values on this cell
198 processor(cell, local_values, values);
199 else
200 {
201 Assert(local_values.size() ==
202 cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell(),
203 ExcMessage("Incorrect size of local_values vector."));
204
205 FullMatrix<double> interpolation(
206 cell.get_fe().n_dofs_per_cell(),
207 cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell());
208
209 cell.get_fe().get_interpolation_matrix(
210 cell.get_dof_handler().get_fe(fe_index), interpolation);
211
212 // do the interpolation to the target space. for historical
213 // reasons, matrices are set to size 0x0 internally even if
214 // we reinit as 4x0, so we have to treat this case specially
215 Vector<number> tmp(cell.get_fe().n_dofs_per_cell());
216 if ((tmp.size() > 0) && (local_values.size() > 0))
217 interpolation.vmult(tmp, local_values);
218
219 // now set the dof values in the global vector
220 processor(cell, tmp, values);
221 }
222 }
223 else
224 // otherwise distribute them to the children
225 {
226 Assert((cell.get_dof_handler().has_hp_capabilities() == false) ||
227 (fe_index != numbers::invalid_fe_index),
229 "You cannot call this function on non-active cells "
230 "of DoFHandler objects unless you provide an explicit "
231 "finite element index because they do not have naturally "
232 "associated finite element spaces associated: degrees "
233 "of freedom are only distributed on active cells for which "
234 "the active FE index has been set."));
235
237 cell.get_dof_handler().get_fe(fe_index);
238 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
239
240 Assert(local_values.size() == dofs_per_cell,
242 ExcVectorDoesNotMatch()));
243 Assert(values.size() == cell.get_dof_handler().n_dofs(),
245 ExcVectorDoesNotMatch()));
246
247 Vector<number> tmp(dofs_per_cell);
248
249 for (unsigned int child = 0; child < cell.n_children(); ++child)
250 {
251 if (tmp.size() > 0)
252 fe.get_prolongation_matrix(child, cell.refinement_case())
253 .vmult(tmp, local_values);
255 *cell.child(child), tmp, values, fe_index, processor);
256 }
257 }
258 }
259
260} // namespace internal
261
262
263
264template <int dim, int spacedim, bool lda>
265template <class OutputVector, typename number>
266void
268 const Vector<number> &local_values,
269 OutputVector &values,
270 const types::fe_index fe_index_,
271 const bool perform_check) const
272{
273 internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
274 *this,
275 local_values,
276 values,
277 fe_index_,
278 [perform_check](const DoFCellAccessor<dim, spacedim, lda> &cell,
279 const Vector<number> &local_values,
280 OutputVector &values) {
281 internal::set_dof_values(cell, local_values, values, perform_check);
282 });
283}
284
285
286template <int dim, int spacedim, bool lda>
287template <class OutputVector, typename number>
288void
291 const Vector<number> &local_values,
292 OutputVector &values,
293 const types::fe_index fe_index_) const
294{
295 internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
296 *this,
297 local_values,
298 values,
299 fe_index_,
301 const Vector<number> &local_values,
302 OutputVector &values) {
303 std::vector<types::global_dof_index> dof_indices(
304 cell.get_fe().n_dofs_per_cell());
305 cell.get_dof_indices(dof_indices);
307 dof_indices,
308 values);
309 });
310}
311
312
313// --------------------------------------------------------------------------
314// explicit instantiations
315#include "dof_accessor_set.inst"
316
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:534
static ::ExceptionBase & ExcNonMatchingElementsSetDofValuesByInterpolation(Number arg1, Number arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
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:243
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:59