Reference documentation for deal.II version GIT d6cf33b98c 2023-09-22 19:50:02+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\}}\)
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 
26 namespace parallel
27 {
28  namespace distributed
29  {
30  namespace experimental
31  {
32  template <int dim, typename VectorType, int spacedim>
34  const DoFHandler<dim, spacedim> &dof)
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(),
50  &dof_handler);
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,
60  ExcMessage(
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(),
77  ExcMessage(
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 (const auto &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 Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void compute_inverse_svd(const double threshold=0.)
Definition: vector.h:110
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int global_dof_index
Definition: types.h:82
unsigned short int fe_index
Definition: types.h:60