Reference documentation for deal.II version GIT c545eda070 2023-01-27 00:25: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 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  cell_data_transfer = std::make_unique<
39  dynamic_cast<
41  const_cast<::Triangulation<dim, spacedim> &>(
43  }
44 
45 
46 
47  template <int dim, typename VectorType, int spacedim>
48  void
51  const VectorType & in,
52  const unsigned int fe_nothing_index)
53  {
54  const unsigned int dofs_per_cell =
55  dof_handler.get_fe_collection().max_dofs_per_cell();
56 
57  Vector<Number> cell_solution(dofs_per_cell);
58  Vector<Number> dummy_cell_solution(dofs_per_cell);
59 
60  for (auto &val : dummy_cell_solution)
61  {
62  val = std::numeric_limits<Number>::infinity();
63  }
64 
65  in.update_ghost_values();
66 
67  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
68  for (const auto &cell : dof_handler.active_cell_iterators())
69  {
70  if ((cell->is_locally_owned()) &&
71  (cell->active_fe_index() != fe_nothing_index))
72  {
73  cell->get_dof_values(in, cell_solution);
74  data_to_transfer.push_back(cell_solution);
75  }
76  else
77  {
78  data_to_transfer.push_back(dummy_cell_solution);
79  }
80  }
81 
82  cell_data_transfer->prepare_for_coarsening_and_refinement(
83  data_to_transfer);
84  }
85 
86 
87 
88  template <int dim, typename VectorType, int spacedim>
89  void
91  const Number & new_value,
92  const AffineConstraints<Number> &affine_constraints,
93  VectorType & out)
94  {
95  const unsigned int dofs_per_cell =
96  dof_handler.get_fe_collection().max_dofs_per_cell();
97  std::vector<Vector<Number>> transferred_data(
98  dof_handler.get_triangulation().n_active_cells(),
99  Vector<Number>(dofs_per_cell));
100 
101  cell_data_transfer->unpack(transferred_data);
102 
103  // Free the memory allocated by data_to_transfer
104  data_to_transfer.clear();
105 
106  for (unsigned int i = 0; i < out.locally_owned_size(); ++i)
107  out.local_element(i) = std::numeric_limits<Number>::infinity();
108 
109  unsigned int cell_i = 0;
110  for (auto const &cell : dof_handler.active_cell_iterators())
111  {
112  if ((cell->is_locally_owned()) &&
113  (transferred_data[cell_i][0] !=
114  std::numeric_limits<Number>::infinity()))
115  {
116  cell->set_dof_values(transferred_data[cell_i], out);
117  }
118  ++cell_i;
119  }
120 
121 
122  // Communicate the results.
123  out.compress(::VectorOperation::insert);
124 
125  // Treat hanging nodes
126  std::vector<types::global_dof_index> dof_indices;
127  std::vector<types::global_dof_index> dofs_map;
128  std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
129  constraint_lines;
130  std::vector<Number> constraint_values;
131  IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
132  for (auto constrained_dof : locally_owned_dofs)
133  if (affine_constraints.is_constrained(constrained_dof))
134  {
135  auto *constraint =
136  affine_constraints.get_constraint_entries(constrained_dof);
137  const unsigned int line_size = constraint->size();
138  bool add_line = false;
139  for (unsigned int i = 0; i < line_size; ++i)
140  {
141  types::global_dof_index constraining_dof =
142  (*constraint)[i].first;
143  // If one of the constraining value is infinity, we need to
144  // reverse the relationship
145  if (out[constraining_dof] ==
146  std::numeric_limits<Number>::infinity())
147  {
148  add_line = true;
149  break;
150  }
151  }
152 
153  if (add_line)
154  {
155  std::vector<std::pair<types::global_dof_index, Number>> line;
156  Number val = out[constrained_dof];
157 
158  for (unsigned int i = 0; i < line_size; ++i)
159  {
160  types::global_dof_index constraining_dof =
161  (*constraint)[i].first;
162 
163  if (out[constraining_dof] ==
164  std::numeric_limits<Number>::infinity())
165  {
166  auto constraining_dof_map_it =
167  std::find(dofs_map.begin(),
168  dofs_map.end(),
169  constraining_dof);
170  if (constraining_dof_map_it == dofs_map.end())
171  {
172  dofs_map.push_back(constraining_dof);
173  }
174  line.push_back((*constraint)[i]);
175  }
176  else
177  {
178  val -=
179  out[constraining_dof] * (*constraint)[i].second;
180  }
181  }
182  constraint_lines.push_back(line);
183  constraint_values.push_back(val);
184  }
185  }
186 
187  // Build a constraint matrix that we invert
188  const unsigned int n_rows = constraint_lines.size();
189  if (n_rows > 0)
190  {
191  const unsigned int n_cols = dofs_map.size();
192  ::LAPACKFullMatrix<Number> constraints_matrix(n_rows, n_cols);
193  ::Vector<Number> constraint_values_vec(n_rows);
194 
195  for (unsigned int i = 0; i < n_rows; ++i)
196  {
197  for (unsigned int j = 0; j < n_cols; ++j)
198  {
199  if (j < constraint_lines[i].size())
200  {
201  auto constraint_it =
202  std::find(dofs_map.begin(),
203  dofs_map.end(),
204  constraint_lines[i][j].first);
205  constraints_matrix(i,
206  constraint_it - dofs_map.begin()) =
207  constraint_lines[i][j].second;
208  }
209  }
210 
211  constraint_values_vec[i] = constraint_values[i];
212  }
213 
214  constraints_matrix.compute_inverse_svd();
215 
216  ::Vector<Number> new_constrained_values(n_cols);
217  constraints_matrix.vmult(new_constrained_values,
218  constraint_values_vec);
219 
220  for (unsigned int i = 0; i < n_cols; ++i)
221  {
222  out[dofs_map[i]] = new_constrained_values[i];
223  }
224  }
225 
226  // Set the value to the newly create DoFs.
227  std::for_each(out.begin(), out.end(), [&](Number &val) {
228  if (val == std::numeric_limits<Number>::infinity())
229  {
230  val = new_value;
231  }
232  });
233  }
234  } // namespace experimental
235  } // namespace distributed
236 } // namespace parallel
237 
238 // explicit instantiations
239 # include "field_transfer.inst"
240 
242 
243 #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
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:109
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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
unsigned int global_dof_index
Definition: types.h:81