Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40: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\}}\)
constrained_linear_operator.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2022 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 
16 #ifndef dealii_constrained_linear_operator_h
17 #define dealii_constrained_linear_operator_h
18 
19 #include <deal.II/base/config.h>
20 
24 
25 
27 
28 
63 template <typename Range, typename Domain, typename Payload>
68 {
69  LinearOperator<Range, Domain, Payload> return_op = exemplar;
70 
71  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
73  ::ExcMessage("The domain and range vectors must be different "
74  "storage locations"));
75 
76  // First, add vector u to v unconditionally and clean up constrained
77  // degrees of freedom later.
78  v += u;
79 
80  const auto &locally_owned_elements = v.locally_owned_elements();
81  for (const auto &line : constraints.get_lines())
82  {
83  const auto i = line.index;
84  if (locally_owned_elements.is_element(i))
85  {
86  v(i) -= u(i);
87  const auto &entries = line.entries;
88  for (types::global_dof_index j = 0; j < entries.size(); ++j)
89  {
90  const auto pos = entries[j].first;
91  v(i) += u(pos) * entries[j].second;
92  }
93  }
94  }
95 
96  v.compress(VectorOperation::add);
97  };
98 
99  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
100  Assert(!::PointerComparison::equal(&v, &u),
101  ::ExcMessage("The domain and range vectors must be different "
102  "storage locations"));
103 
104  // First, add vector u to v unconditionally and clean up constrained
105  // degrees of freedom later.
106  v += u;
107 
108  const auto &locally_owned_elements = v.locally_owned_elements();
109  for (const auto &line : constraints.get_lines())
110  {
111  const auto i = line.index;
112 
113  if (locally_owned_elements.is_element(i))
114  {
115  v(i) -= u(i);
116  }
117 
118  const auto &entries = line.entries;
119  for (types::global_dof_index j = 0; j < entries.size(); ++j)
120  {
121  const auto pos = entries[j].first;
122  if (locally_owned_elements.is_element(pos))
123  v(pos) += u(i) * entries[j].second;
124  }
125  }
126 
127  v.compress(VectorOperation::add);
128  };
129 
130  return_op.vmult = [vmult_add = return_op.vmult_add](Range &v,
131  const Domain &u) {
132  v = 0.;
133  vmult_add(v, u);
134  };
135 
136  return_op.Tvmult = [Tvmult_add = return_op.Tvmult_add](Domain &v,
137  const Range &u) {
138  v = 0.;
139  Tvmult_add(v, u);
140  };
141 
142  return return_op;
143 }
144 
145 
156 template <typename Range, typename Domain, typename Payload>
161 {
162  LinearOperator<Range, Domain, Payload> return_op = exemplar;
163 
164  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
165  const auto &locally_owned_elements = v.locally_owned_elements();
166  for (const auto &line : constraints.get_lines())
167  {
168  const auto i = line.index;
169  if (locally_owned_elements.is_element(i))
170  {
171  v(i) += u(i);
172  }
173  }
174 
175  v.compress(VectorOperation::add);
176  };
177 
178  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
179  const auto &locally_owned_elements = v.locally_owned_elements();
180  for (const auto &line : constraints.get_lines())
181  {
182  const auto i = line.index;
183  if (locally_owned_elements.is_element(i))
184  {
185  v(i) += u(i);
186  }
187  }
188 
189  v.compress(VectorOperation::add);
190  };
191 
192  return_op.vmult = [vmult_add = return_op.vmult_add](Range &v,
193  const Domain &u) {
194  v = 0.;
195  vmult_add(v, u);
196  };
197 
198  return_op.Tvmult = [Tvmult_add = return_op.Tvmult_add](Domain &v,
199  const Range &u) {
200  v = 0.;
201  Tvmult_add(v, u);
202  };
203 
204  return return_op;
205 }
206 
207 
244 template <typename Range, typename Domain, typename Payload>
249 {
250  const auto C = distribute_constraints_linear_operator(constraints, linop);
251  const auto Ct = transpose_operator(C);
252  const auto Id_c = project_to_constrained_linear_operator(constraints, linop);
253  return Ct * linop * C + Id_c;
254 }
255 
256 
290 template <typename Range, typename Domain, typename Payload>
295  const Range &right_hand_side)
296 {
297  PackagedOperation<Range> return_comp;
298 
299  return_comp.reinit_vector = linop.reinit_range_vector;
300 
301  return_comp.apply_add = [&constraints, &linop, &right_hand_side](Range &v) {
302  const auto C = distribute_constraints_linear_operator(constraints, linop);
303  const auto Ct = transpose_operator(C);
304 
305  GrowingVectorMemory<Domain> vector_memory;
306  typename VectorMemory<Domain>::Pointer k(vector_memory);
307  linop.reinit_domain_vector(*k, /*bool fast=*/false);
308  constraints.distribute(*k);
309 
310  v += Ct * (right_hand_side - linop * *k);
311  };
312 
313  return_comp.apply = [apply_add = return_comp.apply_add](Range &v) {
314  v = 0.;
315  apply_add(v);
316  };
317 
318  return return_comp;
319 }
320 
324 
325 #endif
LineRange get_lines() const
void distribute(VectorType &vec) const
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_vector
std::function< void(Range &v)> apply
std::function< void(Range &v)> apply_add
#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:1631
static ::ExceptionBase & ExcMessage(std::string arg1)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
PackagedOperation< Range > constrained_right_hand_side(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop, const Range &right_hand_side)
LinearOperator< Range, Domain, Payload > constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop)
LinearOperator< Range, Domain, Payload > project_to_constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &exemplar)
LinearOperator< Range, Domain, Payload > distribute_constraints_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &exemplar)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
unsigned int global_dof_index
Definition: types.h:82
static bool equal(const T *p1, const T *p2)