Reference documentation for deal.II version 9.2.0
\(\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 - 2019 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 
33 
34 
64 template <typename Range, typename Domain, typename Payload>
69 {
70  LinearOperator<Range, Domain, Payload> return_op = exemplar;
71 
72  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
74  ::ExcMessage("The domain and range vectors must be different "
75  "storage locations"));
76 
77  // First, add vector u to v unconditionally and clean up constrained
78  // degrees of freedom later.
79  v += u;
80 
81  const auto &locally_owned_elements = v.locally_owned_elements();
82  for (const auto &line : constraints.get_lines())
83  {
84  const auto i = line.index;
85  if (locally_owned_elements.is_element(i))
86  {
87  v(i) -= u(i);
88  const auto &entries = line.entries;
89  for (types::global_dof_index j = 0; j < entries.size(); ++j)
90  {
91  const auto pos = entries[j].first;
92  v(i) += u(pos) * entries[j].second;
93  }
94  }
95  }
96 
97  v.compress(VectorOperation::add);
98  };
99 
100  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
101  Assert(!::PointerComparison::equal(&v, &u),
102  ::ExcMessage("The domain and range vectors must be different "
103  "storage locations"));
104 
105  // First, add vector u to v unconditionally and clean up constrained
106  // degrees of freedom later.
107  v += u;
108 
109  const auto &locally_owned_elements = v.locally_owned_elements();
110  for (const auto &line : constraints.get_lines())
111  {
112  const auto i = line.index;
113 
114  if (locally_owned_elements.is_element(i))
115  {
116  v(i) -= u(i);
117  }
118 
119  const auto &entries = line.entries;
120  for (types::global_dof_index j = 0; j < entries.size(); ++j)
121  {
122  const auto pos = entries[j].first;
123  if (locally_owned_elements.is_element(pos))
124  v(pos) += u(i) * entries[j].second;
125  }
126  }
127 
128  v.compress(VectorOperation::add);
129  };
130 
131  // lambda capture expressions are a C++14 feature...
132  const auto vmult_add = return_op.vmult_add;
133  return_op.vmult = [vmult_add](Range &v, const Domain &u) {
134  v = 0.;
135  vmult_add(v, u);
136  };
137 
138  // lambda capture expressions are a C++14 feature...
139  const auto Tvmult_add = return_op.Tvmult_add;
140  return_op.Tvmult = [Tvmult_add](Domain &v, const Range &u) {
141  v = 0.;
142  Tvmult_add(v, u);
143  };
144 
145  return return_op;
146 }
147 
148 
160 template <typename Range, typename Domain, typename Payload>
165 {
166  LinearOperator<Range, Domain, Payload> return_op = exemplar;
167 
168  return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
169  const auto &locally_owned_elements = v.locally_owned_elements();
170  for (const auto &line : constraints.get_lines())
171  {
172  const auto i = line.index;
173  if (locally_owned_elements.is_element(i))
174  {
175  v(i) += u(i);
176  }
177  }
178 
179  v.compress(VectorOperation::add);
180  };
181 
182  return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
183  const auto &locally_owned_elements = v.locally_owned_elements();
184  for (const auto &line : constraints.get_lines())
185  {
186  const auto i = line.index;
187  if (locally_owned_elements.is_element(i))
188  {
189  v(i) += u(i);
190  }
191  }
192 
193  v.compress(VectorOperation::add);
194  };
195 
196  // lambda capture expressions are a C++14 feature...
197  const auto vmult_add = return_op.vmult_add;
198  return_op.vmult = [vmult_add](Range &v, const Domain &u) {
199  v = 0.;
200  vmult_add(v, u);
201  };
202 
203  // lambda capture expressions are a C++14 feature...
204  const auto Tvmult_add = return_op.Tvmult_add;
205  return_op.Tvmult = [Tvmult_add](Domain &v, const Range &u) {
206  v = 0.;
207  Tvmult_add(v, u);
208  };
209 
210  return return_op;
211 }
212 
213 
251 template <typename Range, typename Domain, typename Payload>
256 {
257  const auto C = distribute_constraints_linear_operator(constraints, linop);
258  const auto Ct = transpose_operator(C);
259  const auto Id_c = project_to_constrained_linear_operator(constraints, linop);
260  return Ct * linop * C + Id_c;
261 }
262 
263 
298 template <typename Range, typename Domain, typename Payload>
303  const Range & right_hand_side)
304 {
305  PackagedOperation<Range> return_comp;
306 
307  return_comp.reinit_vector = linop.reinit_range_vector;
308 
309  return_comp.apply_add = [&constraints, &linop, &right_hand_side](Range &v) {
310  const auto C = distribute_constraints_linear_operator(constraints, linop);
311  const auto Ct = transpose_operator(C);
312 
313  GrowingVectorMemory<Domain> vector_memory;
314  typename VectorMemory<Domain>::Pointer k(vector_memory);
315  linop.reinit_domain_vector(*k, /*bool fast=*/false);
316  constraints.distribute(*k);
317 
318  v += Ct * (right_hand_side - linop * *k);
319  };
320 
321  // lambda capture expressions are a C++14 feature...
322  const auto apply_add = return_comp.apply_add;
323  return_comp.apply = [apply_add](Range &v) {
324  v = 0.;
325  apply_add(v);
326  };
327 
328  return return_comp;
329 }
330 
332 
334 
335 #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)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
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:1419
std::function< void(Range &v, const Domain &u)> vmult
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
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)
std::function< void(Range &v, const Domain &u)> vmult_add
void distribute(VectorType &vec) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
static bool equal(const T *p1, const T *p2)