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
21 #include <deal.II/lac/affine_constraints.h>
22 #include <deal.II/lac/linear_operator.h>
23 #include <deal.II/lac/packaged_operation.h>
24
25
26 DEAL_II_NAMESPACE_OPEN
27
28
29 /**
30 * @name Indirectly applying constraints to LinearOperator
31 */
32 //@{
33
34
35 /**
36 * This function takes an AffineConstraints object @p constraints and
37 * an operator exemplar @p exemplar (this exemplar is usually a linear
38 * operator that describes the system matrix - it is only used to create
39 * domain and range vectors of appropriate sizes, its action <tt>vmult</tt>
40 * is never used). A LinearOperator object associated with the "homogeneous
41 * action" of the underlying AffineConstraints object is returned:
42 *
43 * Applying the LinearOperator object on a vector <code>u</code> results in a
44 * vector <code>v</code> that stores the result of calling
45 * AffineConstraints::distribute() on <code>u</code> - with one important
46 * difference: inhomogeneities are not applied, but always treated as 0
47 * instead.
48 *
49 * The LinearOperator object created by this function is primarily used
50 * internally in constrained_linear_operator() to build up a modified system
51 * of linear equations. How to solve a linear system of equations with this
52 * approach is explained in detail in the
53 * @ref constraints
54 * module.
55 *
56 *
57 * @note Currently, this function may not work correctly for distributed data
58 * structures.
59 *
60 * @relatesalso LinearOperator
61 * @ingroup constraints
62 */
63 template <typename Range, typename Domain, typename Payload>
64 LinearOperator<Range, Domain, Payload>
distribute_constraints_linear_operator(const AffineConstraints<typename Range::value_type> & constraints,const LinearOperator<Range,Domain,Payload> & exemplar)65 distribute_constraints_linear_operator(
66 const AffineConstraints<typename Range::value_type> &constraints,
67 const LinearOperator<Range, Domain, Payload> & exemplar)
68 {
69 LinearOperator<Range, Domain, Payload> return_op = exemplar;
70
71 return_op.vmult_add = [&constraints](Range &v, const Domain &u) {
72 Assert(!dealii::PointerComparison::equal(&v, &u),
73 dealii::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(!dealii::PointerComparison::equal(&v, &u),
101 dealii::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
146 /**
147 * Given a AffineConstraints @p constraints and an operator exemplar @p
148 * exemplar, return a LinearOperator that is the projection to the subspace of
149 * constrained degrees of freedom, i.e. all entries of the result vector that
150 * correspond to unconstrained degrees of freedom are set to zero.
151 *
152 *
153 * @relatesalso LinearOperator
154 * @ingroup constraints
155 */
156 template <typename Range, typename Domain, typename Payload>
157 LinearOperator<Range, Domain, Payload>
project_to_constrained_linear_operator(const AffineConstraints<typename Range::value_type> & constraints,const LinearOperator<Range,Domain,Payload> & exemplar)158 project_to_constrained_linear_operator(
159 const AffineConstraints<typename Range::value_type> &constraints,
160 const LinearOperator<Range, Domain, Payload> & exemplar)
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
208 /**
209 * Given a AffineConstraints object @p constraints and a LinearOperator
210 * @p linop, this function creates a LinearOperator object consisting of the
211 * composition of three operations and a regularization:
212 * @code
213 * Ct * linop * C + Id_c;
214 * @endcode
215 * with
216 * @code
217 * C = distribute_constraints_linear_operator(constraints, linop);
218 * Ct = transpose_operator(C);
219 * Id_c = project_to_constrained_linear_operator(constraints, linop);
220 * @endcode
221 * and <code>Id_c</code> is the projection to the subspace consisting of all
222 * vector entries associated with constrained degrees of freedom.
223 *
224 * This LinearOperator object is used together with
225 * constrained_right_hand_side() to build up the following modified system of
226 * linear equations:
227 * @f[
228 * (C^T A C + Id_c) x = C^T (b - A\,k)
229 * @f]
230 * with a given (unconstrained) system matrix $A$, right hand side $b$, and
231 * linear constraints $C$ with inhomogeneities $k$.
232 *
233 * A detailed explanation of this approach is given in the
234 * @ref constraints
235 * module.
236 *
237 *
238 * @note Currently, this function may not work correctly for distributed data
239 * structures.
240 *
241 * @relatesalso LinearOperator
242 * @ingroup constraints
243 */
244 template <typename Range, typename Domain, typename Payload>
245 LinearOperator<Range, Domain, Payload>
constrained_linear_operator(const AffineConstraints<typename Range::value_type> & constraints,const LinearOperator<Range,Domain,Payload> & linop)246 constrained_linear_operator(
247 const AffineConstraints<typename Range::value_type> &constraints,
248 const LinearOperator<Range, Domain, Payload> & linop)
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
257 /**
258 * Given a AffineConstraints object @p constraints, a LinearOperator @p
259 * linop and a right-hand side @p right_hand_side, this function creates a
260 * PackagedOperation that stores the following computation:
261 * @code
262 * Ct * (right_hand_side - linop * k)
263 * @endcode
264 * with
265 * @code
266 * C = distribute_constraints_linear_operator(constraints, linop);
267 * Ct = transpose_operator(C);
268 * @endcode
269 *
270 * This LinearOperator object is used together with
271 * constrained_right_hand_side() to build up the following modified system of
272 * linear equations:
273 * @f[
274 * (C^T A C + Id_c) x = C^T (b - A\,k)
275 * @f]
276 * with a given (unconstrained) system matrix $A$, right hand side $b$, and
277 * linear constraints $C$ with inhomogeneities $k$.
278 *
279 * A detailed explanation of this approach is given in the
280 * @ref constraints
281 * module.
282 *
283 *
284 * @note Currently, this function may not work correctly for distributed data
285 * structures.
286 *
287 * @relatesalso LinearOperator
288 * @ingroup constraints
289 */
290 template <typename Range, typename Domain, typename Payload>
291 PackagedOperation<Range>
constrained_right_hand_side(const AffineConstraints<typename Range::value_type> & constraints,const LinearOperator<Range,Domain,Payload> & linop,const Range & right_hand_side)292 constrained_right_hand_side(
293 const AffineConstraints<typename Range::value_type> &constraints,
294 const LinearOperator<Range, Domain, Payload> & linop,
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
321 //@}
322
323 DEAL_II_NAMESPACE_CLOSE
324
325 #endif
326