Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
constrained_linear_operator.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2018 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/lac/affine_constraints.h>
20 #include <deal.II/lac/linear_operator.h>
21 #include <deal.II/lac/packaged_operation.h>
22 
23 
24 DEAL_II_NAMESPACE_OPEN
25 
26 
31 
32 
62 template <typename Range, typename Domain, typename Payload>
67 {
68  LinearOperator<Range, Domain, Payload> return_op = exemplar;
69 
70  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
72  ::ExcMessage("The domain and range vectors must be different "
73  "storage locations"));
74 
75  // First, add vector u to v unconditionally and clean up constrained
76  // degrees of freedom later.
77  v += u;
78 
79  const auto &locally_owned_elements = v.locally_owned_elements();
80  for (const auto &line : constraints.get_lines())
81  {
82  const auto i = line.index;
83  if (locally_owned_elements.is_element(i))
84  {
85  v(i) -= u(i);
86  const auto &entries = line.entries;
87  for (types::global_dof_index j = 0; j < entries.size(); ++j)
88  {
89  const auto pos = entries[j].first;
90  v(i) += u(pos) * entries[j].second;
91  }
92  }
93  }
94 
95  v.compress(VectorOperation::add);
96  };
97 
98  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
100  ::ExcMessage("The domain and range vectors must be different "
101  "storage locations"));
102 
103  // First, add vector u to v unconditionally and clean up constrained
104  // degrees of freedom later.
105  v += u;
106 
107  const auto &locally_owned_elements = v.locally_owned_elements();
108  for (const auto &line : constraints.get_lines())
109  {
110  const auto i = line.index;
111 
112  if (locally_owned_elements.is_element(i))
113  {
114  v(i) -= u(i);
115  }
116 
117  const auto &entries = line.entries;
118  for (types::global_dof_index j = 0; j < entries.size(); ++j)
119  {
120  const auto pos = entries[j].first;
121  if (locally_owned_elements.is_element(pos))
122  v(pos) += u(i) * entries[j].second;
123  }
124  }
125 
126  v.compress(VectorOperation::add);
127  };
128 
129  // lambda capture expressions are a C++14 feature...
130  const auto vmult_add = return_op.vmult_add;
131  return_op.vmult = [vmult_add](Range &v, const Domain &u) {
132  v = 0.;
133  vmult_add(v, u);
134  };
135 
136  // lambda capture expressions are a C++14 feature...
137  const auto Tvmult_add = return_op.Tvmult_add;
138  return_op.Tvmult = [Tvmult_add](Domain &v, const Range &u) {
139  v = 0.;
140  Tvmult_add(v, u);
141  };
142 
143  return return_op;
144 }
145 
146 
158 template <typename Range, typename Domain, typename Payload>
163 {
164  LinearOperator<Range, Domain, Payload> return_op = exemplar;
165 
166  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
167  const auto &locally_owned_elements = v.locally_owned_elements();
168  for (const auto &line : constraints.get_lines())
169  {
170  const auto i = line.index;
171  if (locally_owned_elements.is_element(i))
172  {
173  v(i) += u(i);
174  }
175  }
176 
177  v.compress(VectorOperation::add);
178  };
179 
180  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
181  const auto &locally_owned_elements = v.locally_owned_elements();
182  for (const auto &line : constraints.get_lines())
183  {
184  const auto i = line.index;
185  if (locally_owned_elements.is_element(i))
186  {
187  v(i) += u(i);
188  }
189  }
190 
191  v.compress(VectorOperation::add);
192  };
193 
194  // lambda capture expressions are a C++14 feature...
195  const auto vmult_add = return_op.vmult_add;
196  return_op.vmult = [vmult_add](Range &v, const Domain &u) {
197  v = 0.;
198  vmult_add(v, u);
199  };
200 
201  // lambda capture expressions are a C++14 feature...
202  const auto Tvmult_add = return_op.Tvmult_add;
203  return_op.Tvmult = [Tvmult_add](Domain &v, const Range &u) {
204  v = 0.;
205  Tvmult_add(v, u);
206  };
207 
208  return return_op;
209 }
210 
211 
249 template <typename Range, typename Domain, typename Payload>
254 {
255  const auto C = distribute_constraints_linear_operator(constraints, linop);
256  const auto Ct = transpose_operator(C);
257  const auto Id_c = project_to_constrained_linear_operator(constraints, linop);
258  return Ct * linop * C + Id_c;
259 }
260 
261 
296 template <typename Range, typename Domain, typename Payload>
301  const Range & right_hand_side)
302 {
303  PackagedOperation<Range> return_comp;
304 
305  return_comp.reinit_vector = linop.reinit_range_vector;
306 
307  return_comp.apply_add = [&constraints, &linop, &right_hand_side](Range &v) {
308  const auto C = distribute_constraints_linear_operator(constraints, linop);
309  const auto Ct = transpose_operator(C);
310 
311  GrowingVectorMemory<Domain> vector_memory;
312  typename VectorMemory<Domain>::Pointer k(vector_memory);
313  linop.reinit_domain_vector(*k, /*bool fast=*/false);
314  constraints.distribute(*k);
315 
316  v += Ct * (right_hand_side - linop * *k);
317  };
318 
319  // lambda capture expressions are a C++14 feature...
320  const auto apply_add = return_comp.apply_add;
321  return_comp.apply = [apply_add](Range &v) {
322  v = 0.;
323  apply_add(v);
324  };
325 
326  return return_comp;
327 }
328 
330 
331 DEAL_II_NAMESPACE_CLOSE
332 
333 #endif
std::function< void(Range &v)> apply
LinearOperator< Range, Domain, Payload > distribute_constraints_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &exemplar)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_vector
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(Domain &v, const Range &u)> Tvmult
const LineRange get_lines() const
LinearOperator< Range, Domain, Payload > constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
std::function< void(Range &v, const Domain &u)> vmult
PackagedOperation< Range > constrained_right_hand_side(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop, const Range &right_hand_side)
std::function< void(Range &v)> apply_add
LinearOperator< Range, Domain, Payload > project_to_constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &exemplar)
unsigned int global_dof_index
Definition: types.h:89
std::function< void(Range &v, const Domain &u)> vmult_add
void distribute(VectorType &vec) const
static bool equal(const T *p1, const T *p2)