Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
63template <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) {
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
156template <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
244template <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
290template <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
const 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
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)
static bool equal(const T *p1, const T *p2)