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
field_transfer.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2021 - 2023 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
17
18#ifdef DEAL_II_WITH_P4EST
19
21
22# include <limits>
23
25
26namespace parallel
27{
28 namespace distributed
29 {
30 namespace experimental
31 {
32 template <int dim, typename VectorType, int spacedim>
35 : dof_handler(dof)
36 {
37 // When coarsening, we want to mimic the behavior of SolutionTransfer
38 // and interpolate from child cells to parent. Define this strategy here
39 // since it is not readily available
40 const auto coarsening_strategy =
41 [this](
42 const typename ::Triangulation<dim, spacedim>::cell_iterator
43 & parent,
44 const std::vector<Vector<Number>> &children_values) {
45 // get the equivalent DoFCellAccessor
46 typename DoFHandler<dim, spacedim>::cell_iterator dof_cell_iterator(
48 parent->level(),
49 parent->index(),
51
52 int fe_index = 0;
54 fe_index = ::internal::hp::DoFHandlerImplementation::
55 dominated_future_fe_on_children<dim, spacedim>(
56 dof_cell_iterator);
57
58 const auto &fe = dof_handler.get_fe(fe_index);
59 Assert(fe.n_dofs_per_cell() > 0,
61 "Cannot coarsen onto a FiniteElement with no DoFs."));
62 AssertDimension(dof_cell_iterator->n_children(),
63 children_values.size());
64
65 const auto child_iterators = dof_cell_iterator->child_iterators();
66 const unsigned int n_children_with_fe_nothing =
67 std::count_if(child_iterators.begin(),
68 child_iterators.end(),
69 [](const auto &child_cell) {
70 return child_cell->get_fe().n_dofs_per_cell() ==
71 0;
72 });
73
74 Assert(
75 n_children_with_fe_nothing == 0 ||
76 n_children_with_fe_nothing == dof_cell_iterator->n_children(),
78 "Coarsening is only supported for parent cells where either all"
79 " or none of the child cells are FE_Nothing."));
80
81 // in case all children are FE_Nothing there is nothing to
82 // interpolate and we just return the first entry from the children
83 // values (containing invalid entries)
84 if (n_children_with_fe_nothing == dof_cell_iterator->n_children())
85 {
86 return children_values[0];
87 }
88
89 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
90 Vector<Number> tmp(dofs_per_cell);
91 Vector<Number> interpolated_values(dofs_per_cell);
92
93 // Otherwise, perform the actual interpolation here. Due to the
94 // assert above, we know that all child cells have data to
95 // interpolate.
96 for (unsigned int child = 0;
97 child < dof_cell_iterator->n_children();
98 ++child)
99 {
100 // interpolate the previously stored values on a child to the
101 // mother cell
102 fe.get_restriction_matrix(child,
103 dof_cell_iterator->refinement_case())
104 .vmult(tmp, children_values[child]);
105
106 // and add up or set them in the output vector
107 for (unsigned int i = 0; i < dofs_per_cell; ++i)
108 if (fe.restriction_is_additive(i))
109 interpolated_values(i) += tmp(i);
110 else if (tmp(i) != Number())
111 interpolated_values(i) = tmp(i);
112 }
113
114 return interpolated_values;
115 };
116
117 cell_data_transfer = std::make_unique<
119 dynamic_cast<
121 const_cast<::Triangulation<dim, spacedim> &>(
123 false,
124 &::AdaptationStrategies::Refinement::
125 preserve<dim, spacedim, Vector<Number>>,
126 coarsening_strategy);
127 }
128
129
130
131 template <int dim, typename VectorType, int spacedim>
132 void
135 const VectorType & in,
136 const unsigned int fe_nothing_index)
137 {
138 const unsigned int dofs_per_cell =
139 dof_handler.get_fe_collection().max_dofs_per_cell();
140
141 Vector<Number> cell_solution(dofs_per_cell);
142 Vector<Number> dummy_cell_solution(dofs_per_cell);
143
144 for (auto &val : dummy_cell_solution)
145 {
146 val = std::numeric_limits<Number>::infinity();
147 }
148
149 in.update_ghost_values();
150
151 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
152 for (const auto &cell : dof_handler.active_cell_iterators())
153 {
154 if ((cell->is_locally_owned()) &&
155 (cell->active_fe_index() != fe_nothing_index))
156 {
157 cell->get_dof_values(in, cell_solution);
158 data_to_transfer.push_back(cell_solution);
159 }
160 else
161 {
162 data_to_transfer.push_back(dummy_cell_solution);
163 }
164 }
165
166 cell_data_transfer->prepare_for_coarsening_and_refinement(
167 data_to_transfer);
168 }
169
170
171
172 template <int dim, typename VectorType, int spacedim>
173 void
175 const Number & new_value,
176 const AffineConstraints<Number> &affine_constraints,
177 VectorType & out)
178 {
179 const unsigned int dofs_per_cell =
180 dof_handler.get_fe_collection().max_dofs_per_cell();
181 std::vector<Vector<Number>> transferred_data(
182 dof_handler.get_triangulation().n_active_cells(),
183 Vector<Number>(dofs_per_cell));
184
185 cell_data_transfer->unpack(transferred_data);
186
187 // Free the memory allocated by data_to_transfer
188 data_to_transfer.clear();
189
190 for (unsigned int i = 0; i < out.locally_owned_size(); ++i)
191 out.local_element(i) = std::numeric_limits<Number>::infinity();
192
193 unsigned int cell_i = 0;
194 for (auto const &cell : dof_handler.active_cell_iterators())
195 {
196 if ((cell->is_locally_owned()) &&
197 (transferred_data[cell_i][0] !=
198 std::numeric_limits<Number>::infinity()))
199 {
200 cell->set_dof_values(transferred_data[cell_i], out);
201 }
202 ++cell_i;
203 }
204
205
206 // Communicate the results.
207 out.compress(VectorOperation::insert);
208
209 // Treat hanging nodes
210 std::vector<types::global_dof_index> dof_indices;
211 std::vector<types::global_dof_index> dofs_map;
212 std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
213 constraint_lines;
214 std::vector<Number> constraint_values;
215 IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
216 for (auto constrained_dof : locally_owned_dofs)
217 if (affine_constraints.is_constrained(constrained_dof))
218 {
219 auto *constraint =
220 affine_constraints.get_constraint_entries(constrained_dof);
221 const unsigned int line_size = constraint->size();
222 bool add_line = false;
223 for (unsigned int i = 0; i < line_size; ++i)
224 {
225 types::global_dof_index constraining_dof =
226 (*constraint)[i].first;
227 // If one of the constraining value is infinity, we need to
228 // reverse the relationship
229 if (out[constraining_dof] ==
230 std::numeric_limits<Number>::infinity())
231 {
232 add_line = true;
233 break;
234 }
235 }
236
237 if (add_line)
238 {
239 std::vector<std::pair<types::global_dof_index, Number>> line;
240 Number val = out[constrained_dof];
241
242 for (unsigned int i = 0; i < line_size; ++i)
243 {
244 types::global_dof_index constraining_dof =
245 (*constraint)[i].first;
246
247 if (out[constraining_dof] ==
248 std::numeric_limits<Number>::infinity())
249 {
250 auto constraining_dof_map_it =
251 std::find(dofs_map.begin(),
252 dofs_map.end(),
253 constraining_dof);
254 if (constraining_dof_map_it == dofs_map.end())
255 {
256 dofs_map.push_back(constraining_dof);
257 }
258 line.push_back((*constraint)[i]);
259 }
260 else
261 {
262 val -=
263 out[constraining_dof] * (*constraint)[i].second;
264 }
265 }
266 constraint_lines.push_back(line);
267 constraint_values.push_back(val);
268 }
269 }
270
271 // Build a constraint matrix that we invert
272 const unsigned int n_rows = constraint_lines.size();
273 if (n_rows > 0)
274 {
275 const unsigned int n_cols = dofs_map.size();
276 ::LAPACKFullMatrix<Number> constraints_matrix(n_rows, n_cols);
277 ::Vector<Number> constraint_values_vec(n_rows);
278
279 for (unsigned int i = 0; i < n_rows; ++i)
280 {
281 for (unsigned int j = 0; j < n_cols; ++j)
282 {
283 if (j < constraint_lines[i].size())
284 {
285 auto constraint_it =
286 std::find(dofs_map.begin(),
287 dofs_map.end(),
288 constraint_lines[i][j].first);
289 constraints_matrix(i,
290 constraint_it - dofs_map.begin()) =
291 constraint_lines[i][j].second;
292 }
293 }
294
295 constraint_values_vec[i] = constraint_values[i];
296 }
297
298 constraints_matrix.compute_inverse_svd();
299
300 ::Vector<Number> new_constrained_values(n_cols);
301 constraints_matrix.vmult(new_constrained_values,
302 constraint_values_vec);
303
304 for (unsigned int i = 0; i < n_cols; ++i)
305 {
306 out[dofs_map[i]] = new_constrained_values[i];
307 }
308 }
309
310 // Set the value to the newly create DoFs.
311 std::for_each(out.begin(), out.end(), [&](Number &val) {
312 if (val == std::numeric_limits<Number>::infinity())
313 {
314 val = new_value;
315 }
316 });
317 }
318 } // namespace experimental
319 } // namespace distributed
320} // namespace parallel
321
322// explicit instantiations
323# include "field_transfer.inst"
324
326
327#endif
bool is_constrained(const size_type line_n) const
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void compute_inverse_svd(const double threshold=0.)
std::unique_ptr< CellDataTransfer< dim, spacedim, std::vector< Vector< Number > > > > cell_data_transfer
void prepare_for_coarsening_and_refinement(const VectorType &in, const unsigned int fe_nothing_index)
FieldTransfer(const DoFHandler< dim, spacedim > &dof_handler)
const DoFHandler< dim, spacedim > & dof_handler
void interpolate(const Number &new_value, const AffineConstraints< Number > &affine_constraints, VectorType &out)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator