DUNE PDELab (git)

blockoffdiagonalwrapper.hh
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
6#include <dune/pdelab/common/intersectiontype.hh>
7#include <dune/pdelab/localoperator/blockdiagonalwrapper.hh>
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>
15 class ZeroViewWrapper
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:
58 using View = ConstAliasedVectorView<Container, LocalFunctionSpaceCache>;
59 using ElementType = typename View::ElementType;
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.
164 impl::BlockDiagonalAccumulationViewWrapper<MAT> view_ns(mat_ns, true);
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>
170 void jacobian_apply_skeleton(const IG& ig,
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
195 impl::BlockDiagonalAccumulationViewWrapper<Y> view_s_on(y_s, true);
196 impl::BlockDiagonalAccumulationViewWrapper<Y> view_n_off(y_n, false);
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>
203 void jacobian_apply_skeleton(const IG& ig,
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
228 impl::BlockDiagonalAccumulationViewWrapper<Y> view_s_on(y_s, true);
229 impl::BlockDiagonalAccumulationViewWrapper<Y> view_n_off(y_n, false);
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
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
A local operator that accumulates the off block diagonal.
Definition: blockoffdiagonalwrapper.hh:112
BlockOffDiagonalLocalOperatorWrapper(const LocalOperator &localOperator)
Construct new instance of class.
Definition: blockoffdiagonalwrapper.hh:135
BlockOffDiagonalLocalOperatorWrapper(const BlockOffDiagonalLocalOperatorWrapper &other)
Copy constructor.
Definition: blockoffdiagonalwrapper.hh:140
Default flags for all local operators.
Definition: flags.hh:19
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)