Reference documentation for deal.II version GIT relicensing-769-g16d39143b0 2024-05-28 17:50: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\}}\)
Loading...
Searching...
No Matches
constrained_linear_operator.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_constrained_linear_operator_h
16#define dealii_constrained_linear_operator_h
17
18#include <deal.II/base/config.h>
19
23
24
26
27
62template <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 return_op.vmult = [vmult_add = return_op.vmult_add](Range &v,
130 const Domain &u) {
131 v = 0.;
132 vmult_add(v, u);
133 };
134
135 return_op.Tvmult = [Tvmult_add = return_op.Tvmult_add](Domain &v,
136 const Range &u) {
137 v = 0.;
138 Tvmult_add(v, u);
139 };
140
141 return return_op;
142}
143
144
155template <typename Range, typename Domain, typename Payload>
160{
161 LinearOperator<Range, Domain, Payload> return_op = exemplar;
162
163 return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
164 const auto &locally_owned_elements = v.locally_owned_elements();
165 for (const auto &line : constraints.get_lines())
166 {
167 const auto i = line.index;
168 if (locally_owned_elements.is_element(i))
169 {
170 v(i) += u(i);
171 }
172 }
173
174 v.compress(VectorOperation::add);
175 };
176
177 return_op.Tvmult_add = [&constraints](Domain &v, const Range &u) {
178 const auto &locally_owned_elements = v.locally_owned_elements();
179 for (const auto &line : constraints.get_lines())
180 {
181 const auto i = line.index;
182 if (locally_owned_elements.is_element(i))
183 {
184 v(i) += u(i);
185 }
186 }
187
188 v.compress(VectorOperation::add);
189 };
190
191 return_op.vmult = [vmult_add = return_op.vmult_add](Range &v,
192 const Domain &u) {
193 v = 0.;
194 vmult_add(v, u);
195 };
196
197 return_op.Tvmult = [Tvmult_add = return_op.Tvmult_add](Domain &v,
198 const Range &u) {
199 v = 0.;
200 Tvmult_add(v, u);
201 };
202
203 return return_op;
204}
205
206
243template <typename Range, typename Domain, typename Payload>
248{
249 const auto C = distribute_constraints_linear_operator(constraints, linop);
250 const auto Ct = transpose_operator(C);
251 const auto Id_c = project_to_constrained_linear_operator(constraints, linop);
252 return Ct * linop * C + Id_c;
253}
254
255
289template <typename Range, typename Domain, typename Payload>
294 const Range &right_hand_side)
295{
296 PackagedOperation<Range> return_comp;
297
298 return_comp.reinit_vector = linop.reinit_range_vector;
299
300 return_comp.apply_add = [&constraints, &linop, &right_hand_side](Range &v) {
301 const auto C = distribute_constraints_linear_operator(constraints, linop);
302 const auto Ct = transpose_operator(C);
303
304 GrowingVectorMemory<Domain> vector_memory;
305 typename VectorMemory<Domain>::Pointer k(vector_memory);
306 linop.reinit_domain_vector(*k, /*bool fast=*/false);
307 constraints.distribute(*k);
308
309 v += Ct * (right_hand_side - linop * *k);
310 };
311
312 return_comp.apply = [apply_add = return_comp.apply_add](Range &v) {
313 v = 0.;
314 apply_add(v);
315 };
316
317 return return_comp;
318}
319
323
324#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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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)