dune-pdelab 2.7-git
Loading...
Searching...
No Matches
blockoffdiagonalwrapper.hh
Go to the documentation of this file.
1// -*- tab-width: 2; indent-tabs-mode: nil -*-
2// vi: set et ts=2 sw=2 sts=2:
3#ifndef DUNE_PDELAB_LOCALOPERATOR_BLOCKOFFDIAGONALWRAPPER_HH
4#define DUNE_PDELAB_LOCALOPERATOR_BLOCKOFFDIAGONALWRAPPER_HH
5
8namespace Dune {
9 namespace PDELab {
10
11 namespace impl {
12
13 // This can be used to get a vector view that returns a zero coefficient.
14 template <typename View>
16 {
17 public:
18 using Container = typename View::Container;
19 using ElementType = typename View::value_type;
20 using SizeType = typename View::size_type;
21
22 // If zero is set to false this class will forward all container
23 // accesses to the vector view that is passed as first argument. This
24 // means it will basically behave the same way as the view itself.
25 //
26 // If you set zero to false it will return 0.0 when it is asked for a
27 // coefficient.
28 //
29 // Use case: The methods in the localoperator interface sometimes get
30 // multiple coefficient views that need to have the same type (e.g. x_s
31 // and x_n for the ansatz function in skeleton terms). This can be used
32 // to 'null' one of those vectors without actually changing any values
33 // in memory.
34 ZeroViewWrapper(const View& view, bool zero)
35 : _view(view), _zero(zero), _zeroCoefficient(0.0)
36 {}
37
38 template <typename LFS>
39 const ElementType& operator()(const LFS& lfs, SizeType i) const
40 {
41 if (_zero)
42 return _zeroCoefficient;
43 else
44 return _view.container()(lfs, i);
45 }
46
47 private:
48 const View& _view;
49 bool _zero;
50 ElementType _zeroCoefficient;
51 };
52
53 // Interfaces look different in the fast-DG case
54 template <typename Container, typename LocalFunctionSpaceCache>
55 class ZeroViewWrapper<AliasedVectorView<Container, LocalFunctionSpaceCache>>
56 {
57 public:
60 using SizeType = typename View::size_type;
61
62 ZeroViewWrapper(const View& view, bool zero)
63 : _view(view), _zero(zero), _zeroCoefficient(0.0)
64 {}
65
66 template <typename LFS>
67 const ElementType& operator()(const LFS& lfs, SizeType i) const
68 {
69 if (_zero)
70 return _zeroCoefficient;
71 else
72 return _view(lfs, i);
73 }
74
75 const ElementType* data() const
76 {
77 // Note: There is no principal problem in implementing this. Create
78 // an array of ElementType that has the correct size and contains
79 // only zeros. This was not implemted since there was no way of
80 // testing the implementation. Better to have a clear error message
81 // than a delicate implementation bug.
82 DUNE_THROW(Dune::Exception, "So far the ZeroViewWrapper does not support fast DG local operators using the data() method to access coefficients. .");
83 }
84
85 private:
86 const View& _view;
87 bool _zero;
88 ElementType _zeroCoefficient;
89 };
90
91
92 } // namespace impl
93
109 template <typename LocalOperator>
112 {
113 public:
114 // We only want to assemble the off-diagonal blocks so we only need
115 // skeleton pattern
116 static constexpr bool doPatternSkeleton = true;
117
118 // For assembling the off-diagonal blocks we only need alpha skeleton
119 static constexpr bool doAlphaSkeleton = LocalOperator::doAlphaSkeleton;
120
121 // If the underlying lop is linear, this one will be linear too
122 static constexpr bool isLinear = LocalOperator::isLinear;
123
124 // This wrapper is designed to use two sided assembly. If it was just
125 // about assembling the whole diagonal block matrix one sided assembly
126 // would be more efficient. For the implementation of matrix-free
127 // preconditioner we need to assemble a diagonal block locally for a
128 // given cell so we need to assemble in a two sided fashion.
129 static constexpr bool doSkeletonTwoSided = true;
130
135 BlockOffDiagonalLocalOperatorWrapper(const LocalOperator& localOperator)
136 : _localOperator(localOperator)
137 {}
138
141 : _localOperator(other._localOperator)
142 {}
143
144 // define sparsity pattern connecting self and neighbor dofs
145 template<typename LFSU, typename LFSV, typename LocalPattern>
146 void pattern_skeleton (const LFSU& lfsu_s, const LFSV& lfsv_s, const LFSU& lfsu_n, const LFSV& lfsv_n,
147 LocalPattern& pattern_sn,
148 LocalPattern& pattern_ns) const
149 {
150 _localOperator.pattern_skeleton (lfsu_s, lfsv_s, lfsu_n, lfsv_n, pattern_sn, pattern_ns);
151 }
152
153 template<typename IG, typename LFSU, typename X, typename LFSV, typename MAT>
154 void jacobian_skeleton (const IG& ig,
155 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
156 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
157 MAT& mat_ss, MAT& mat_sn,
158 MAT& mat_ns, MAT& mat_nn) const
159 {
160 // Since we do two sided assembly (visiting each intersection twice) we
161 // have options. We could assemble either mat_sn or mat_ns. Both lead
162 // to the same solution. In order to be consistent with the choice for
163 // jacobian_apply_skeleton we will assemble m_ns.
165 impl::BlockDiagonalAccumulationViewWrapper<MAT> view_other(mat_ns, false);
166 _localOperator.jacobian_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, view_other, view_other, view_ns, view_other);
167 }
168
169 template<typename IG, typename LFSU, typename Z, typename LFSV, typename Y>
171 const LFSU& lfsu_s, const Z& z_s, const LFSV& lfsv_s,
172 const LFSU& lfsu_n, const Z& z_n, const LFSV& lfsv_n,
173 Y& y_s, Y& y_n) const
174 {
175 // This is more tricky. For a full Jacobian apply you would do:
176 //
177 // y_s = A_ss z_s + A_ns z_n
178 // y_n = A_sn z_s + A_nn z_n
179 //
180 // Instead we want:
181 //
182 // (1) y_s = A_ns z_n
183 // (2) y_n = A_sn z_s
184 //
185 // Since we do two sided assembly we only need to assemble one of these
186 // two. For the matrix free preconditioner we need to apply the
187 // jacobian locally for one block so we need to implement equation (1)
188 // to get all contributions of other cell on the current one.
189
190 // Set input coefficients z_s to zero
191 impl::ZeroViewWrapper<Z> z_zero(z_s, true);
192 impl::ZeroViewWrapper<Z> z_neigh(z_n, false);
193
194 // Only accumulate in y_s
197
198 // Apply Jacobian
199 Dune::PDELab::impl::jacobianApplySkeleton(_localOperator, ig, lfsu_s, z_zero, lfsv_s, lfsu_n, z_neigh, lfsv_n, view_s_on, view_n_off);
200 }
201
202 template<typename IG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
204 const LFSU& lfsu_s, const X& x_s, const Z& z_s, const LFSV& lfsv_s,
205 const LFSU& lfsu_n, const X& x_n, const Z& z_n, const LFSV& lfsv_n,
206 Y& y_s, Y& y_n) const
207 {
208 // This is more tricky. For a full Jacobian apply you would do:
209 //
210 // y_s = A_ss z_s + A_ns z_n
211 // y_n = A_sn z_s + A_nn z_n
212 //
213 // Instead we want:
214 //
215 // (1) y_s = A_ns z_n
216 // (2) y_n = A_sn z_s
217 //
218 // Since we do two sided assembly we only need to assemble one of these
219 // two. For the matrix free preconditioner we need to apply the
220 // jacobian locally for one block so we need to implement equation (1)
221 // to get all contributions of other cell on the current one.
222
223 // Set input coefficients z_s to zero
224 impl::ZeroViewWrapper<Z> z_zero(z_s, true);
225 impl::ZeroViewWrapper<Z> z_neigh(z_n, false);
226
227 // Only accumulate in y_s
230
231 // Apply Jacobian
232 Dune::PDELab::impl::jacobianApplySkeleton(_localOperator, ig, lfsu_s, x_s, z_zero, lfsv_s, lfsu_n, x_n, z_neigh, lfsv_n, view_s_on, view_n_off);
233 }
234
235 private:
237 const LocalOperator& _localOperator;
238 };
239
240 } // namespace PDELab
241} // namespace Dune
242
243#endif
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
std::enable_if_t< LOP::isLinear > jacobianApplySkeleton(const LOP &lop, const IG &ig, const LFSU &lfsu_s, const X &z_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &z_n, const LFSV &lfsv_n, Y &y_s, Y &y_n)
Definition: jacobianapplyhelper.hh:64
Definition: aliasedvectorview.hh:18
Container::size_type size_type
Definition: aliasedvectorview.hh:24
Container::E ElementType
Definition: aliasedvectorview.hh:23
Definition: aliasedvectorview.hh:128
Definition: blockdiagonalwrapper.hh:19
Definition: blockoffdiagonalwrapper.hh:16
ZeroViewWrapper(const View &view, bool zero)
Definition: blockoffdiagonalwrapper.hh:34
typename View::Container Container
Definition: blockoffdiagonalwrapper.hh:18
typename View::value_type ElementType
Definition: blockoffdiagonalwrapper.hh:19
const ElementType & operator()(const LFS &lfs, SizeType i) const
Definition: blockoffdiagonalwrapper.hh:39
typename View::size_type SizeType
Definition: blockoffdiagonalwrapper.hh:20
typename View::size_type SizeType
Definition: blockoffdiagonalwrapper.hh:60
ZeroViewWrapper(const View &view, bool zero)
Definition: blockoffdiagonalwrapper.hh:62
typename View::ElementType ElementType
Definition: blockoffdiagonalwrapper.hh:59
const ElementType & operator()(const LFS &lfs, SizeType i) const
Definition: blockoffdiagonalwrapper.hh:67
const ElementType * data() const
Definition: blockoffdiagonalwrapper.hh:75
A local operator that accumulates the off block diagonal.
Definition: blockoffdiagonalwrapper.hh:112
static constexpr bool doSkeletonTwoSided
Definition: blockoffdiagonalwrapper.hh:129
static constexpr bool isLinear
Definition: blockoffdiagonalwrapper.hh:122
static constexpr bool doAlphaSkeleton
Definition: blockoffdiagonalwrapper.hh:119
BlockOffDiagonalLocalOperatorWrapper(const LocalOperator &localOperator)
Construct new instance of class.
Definition: blockoffdiagonalwrapper.hh:135
void jacobian_apply_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const Z &z_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const Z &z_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
Definition: blockoffdiagonalwrapper.hh:203
static constexpr bool doPatternSkeleton
Definition: blockoffdiagonalwrapper.hh:116
void jacobian_apply_skeleton(const IG &ig, const LFSU &lfsu_s, const Z &z_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const Z &z_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
Definition: blockoffdiagonalwrapper.hh:170
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, MAT &mat_ss, MAT &mat_sn, MAT &mat_ns, MAT &mat_nn) const
Definition: blockoffdiagonalwrapper.hh:154
BlockOffDiagonalLocalOperatorWrapper(const BlockOffDiagonalLocalOperatorWrapper &other)
Copy constructor.
Definition: blockoffdiagonalwrapper.hh:140
void pattern_skeleton(const LFSU &lfsu_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const LFSV &lfsv_n, LocalPattern &pattern_sn, LocalPattern &pattern_ns) const
Definition: blockoffdiagonalwrapper.hh:146
Default flags for all local operators.
Definition: flags.hh:19