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