dune-pdelab 2.7-git
Loading...
Searching...
No Matches
blockdiagonalwrapper.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
4#ifndef DUNE_PDELAB_LOCALOPERATOR_BLOCKDIAGONALWRAPPER_HH
5#define DUNE_PDELAB_LOCALOPERATOR_BLOCKDIAGONALWRAPPER_HH
6
9
10namespace Dune {
11 namespace PDELab {
12
13 namespace impl{
14
15 // This wraps an accumulation view and only accumulates if diagonal is
16 // set to true
17 template <typename AccumulationView>
19 {
20 public:
21
22 BlockDiagonalAccumulationViewWrapper(AccumulationView& view, bool diagonal)
23 : _view(view), _diagonal(diagonal)
24 {}
25
26 const auto weight()
27 {
28 return _view.weight();
29 }
30
31 auto data()
32 {
33 return _view.data();
34 }
35
36 const auto data() const
37 {
38 return _view.data();
39 }
40
41 template <typename LFSU, typename I, typename Value>
42 void accumulate(const LFSU& lfsu, I i, Value value)
43 {
44 if (_diagonal){
45 _view.accumulate(lfsu, i, value);
46 }
47 }
48
49 template <typename LFSU, typename I, typename LFSV, typename J, typename Value>
50 void accumulate(const LFSU& lfsu, I i, const LFSV& lfsv, J j, Value value)
51 {
52 if (_diagonal){
53 _view.accumulate(lfsu, i, lfsv, j, value);
54 }
55 }
56
57 private:
58 AccumulationView& _view;
59 bool _diagonal;
60 };
61 } // namespace impl
62
63
79 template <class LocalOperator>
82 {
83 public:
84 // We only want to assemble the diagonal blocks so we only need volume
85 // pattern.
86 static constexpr bool doPatternVolume = true;
87
88 // For assembling the diagonal block correctly we might need to call
89 // volume, skeleton and boundary
90 static constexpr bool doAlphaVolume = LocalOperator::doAlphaVolume;
91 static constexpr bool doAlphaSkeleton = LocalOperator::doAlphaSkeleton;
92 static constexpr bool doAlphaBoundary = LocalOperator::doAlphaBoundary;
93
94 // If the underlying lop is linear, this one will be linear too
95 static constexpr bool isLinear = LocalOperator::isLinear;
96
97 // This wrapper is designed to use two sided assembly. If it was just
98 // about assembling the whole diagonal block matrix one sided assembly
99 // would be more efficient. For the implementation of matrix-free
100 // preconditioner we need to assemble a diagonal block locally for a
101 // given cell so we need to assemble in a two sided fashion.
102 static constexpr bool doSkeletonTwoSided = true;
103
108 BlockDiagonalLocalOperatorWrapper(const LocalOperator& localOperator)
109 : _localOperator(localOperator)
110 {}
111
114 : _localOperator(other._localOperator)
115 {}
116
117 // define sparsity pattern of operator representation
118 template<typename LFSU, typename LFSV, typename LocalPattern>
119 void pattern_volume (const LFSU& lfsu, const LFSV& lfsv,
120 LocalPattern& pattern) const
121 {
122 _localOperator.pattern_volume(lfsu, lfsv, pattern);
123 }
124
125 template<typename EG, typename LFSU, typename X, typename LFSV, typename MAT>
126 void jacobian_volume (const EG& eg,
127 const LFSU& lfsu,
128 const X& x,
129 const LFSV& lfsv,
130 MAT& mat) const
131 {
132 _localOperator.jacobian_volume(eg,lfsu,x,lfsv,mat);
133 }
134
135 template<typename IG, typename LFSU, typename X, typename LFSV, typename MAT>
136 void jacobian_skeleton (const IG& ig,
137 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
138 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
139 MAT& mat_ss, MAT& mat_sn,
140 MAT& mat_ns, MAT& mat_nn) const
141 {
142 // In the end this jacobian_skeleton only accumulates the self-self
143 // part which corresponds to a diagonal block. The localoperator uses
144 // two sided assembly, so we can discard the neighbor-neighbor block
145 // since it will be accumulated through another cell.
147 impl::BlockDiagonalAccumulationViewWrapper<MAT> view_other(mat_ss, false);
148 _localOperator.jacobian_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, view_ss, view_other, view_other, view_other);
149 }
150
151 template<typename IG, typename LFSU, typename X, typename LFSV, typename MAT>
152 void jacobian_boundary (const IG& ig,
153 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
154 MAT& mat_ss) const
155 {
156 _localOperator.jacobian_boundary(ig, lfsu_s, x_s, lfsv_s, mat_ss);
157 }
158
159
160 template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
161 void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& z, const LFSV& lfsv, Y& y) const
162 {
163 Dune::PDELab::impl::jacobianApplyVolume(_localOperator, eg, lfsu, z, lfsv, y);
164 }
165 template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
166 void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
167 {
168 Dune::PDELab::impl::jacobianApplyVolume(_localOperator, eg, lfsu, x, z, lfsv, y);
169 }
170
171
172 template<typename IG, typename LFSU, typename Z, typename LFSV, typename Y>
174 const LFSU& lfsu_s, const Z& z_s, const LFSV& lfsv_s,
175 const LFSU& lfsu_n, const Z& z_n, const LFSV& lfsv_n,
176 Y& y_s, Y& y_n) const
177 {
180 // Note for the application of only the diagonal block we need to pass
181 // a coefficient vector 0 for the test function into the original
182 // lop. Otherwise it will also apply one of the off-diagonal blocks.
183 Z z_zero(z_n.size(), 0.0);
184 Dune::PDELab::impl::jacobianApplySkeleton(_localOperator, ig, lfsu_s, z_s, lfsv_s, lfsu_n, z_zero, lfsv_n, view_s, view_n);
185 }
186 template<typename IG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
188 const LFSU& lfsu_s, const X& x_s, const Z& z_s, const LFSV& lfsv_s,
189 const LFSU& lfsu_n, const X& x_n, const Z& z_n, const LFSV& lfsv_n,
190 Y& y_s, Y& y_n) const
191 {
194 // Note for the application of only the diagonal block we need to pass
195 // a coefficient vector 0 for the test function into the original
196 // lop. Otherwise it will also apply one of the off-diagonal blocks.
197 Z z_zero(z_n.size(), 0.0);
198 Dune::PDELab::impl::jacobianApplySkeleton(_localOperator, ig, lfsu_s, x_s, z_s, lfsv_s, lfsu_n, x_n, z_zero, lfsv_n, view_s, view_n);
199 }
200
201 template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
203 const LFSU& lfsu_s, const X& z_s, const LFSV& lfsv_s,
204 Y& y_s) const
205 {
206 Dune::PDELab::impl::jacobianApplyBoundary(_localOperator, ig, lfsu_s, z_s, lfsv_s, y_s);
207 }
208
209 template<typename IG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
211 const LFSU& lfsu_s, const X& x_s, const Z& z_s, const LFSV& lfsv_s,
212 Y& y_s) const
213 {
214 Dune::PDELab::impl::jacobianApplyBoundary(_localOperator, ig, lfsu_s, x_s, z_s, lfsv_s, y_s);
215 }
216
217
218 private:
220 const LocalOperator& _localOperator;
221 };
222
223
225 template <typename LocalOperator, typename EG, typename LFSU, typename X, typename LFSV, typename MAT>
226 void assembleLocalDiagonalBlock(const LocalOperator& localOperator,
227 const EG& eg,
228 const LFSU& lfsu,
229 const X& x,
230 const LFSV& lfsv,
231 MAT& mat)
232 {
233 // Assemble the volume part
234 localOperator.jacobian_volume(eg, lfsu, x, lfsv, mat);
235
236 // Iterate over the intersections
237 auto entitySet = lfsu.gridFunctionSpace().entitySet();
238 std::size_t intersectionIndex = 0;
239 for (const auto& is : intersections(lfsu.gridFunctionSpace().gridView(), eg.entity()))
240 {
241 Dune::PDELab::IntersectionGeometry<std::decay_t<decltype(is)>> ig(is, intersectionIndex++);
242 auto intersectionData = classifyIntersection(entitySet, is);
243 auto intersectionType = std::get<0>(intersectionData);
244
245 // Assemble the intersection part
246 switch (intersectionType){
248 localOperator.jacobian_skeleton(ig, lfsu, x, lfsv, lfsu, x, lfsv, mat, mat, mat, mat);
249 break;
251 localOperator.jacobian_boundary(ig, lfsu, x, lfsv, mat);
252 break;
253 default:
254 break;
255 }
256 }
257 }
258
260 template<typename LocalOperator, typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
261 void applyLocalDiagonalBlock(const LocalOperator& localOperator,
262 const EG& eg,
263 const LFSU& lfsu,
264 const X& x,
265 const Z& z,
266 const LFSV& lfsv,
267 Y& y)
268 {
269 // Apply the volume part
270 if (LocalOperator::isLinear)
271 Dune::PDELab::impl::jacobianApplyVolume(localOperator, eg, lfsu, z, lfsv, y);
272 else
273 Dune::PDELab::impl::jacobianApplyVolume(localOperator, eg, lfsu, x, z, lfsv, y);
274
275 // Iterate over the intersections
276 auto entitySet = lfsu.gridFunctionSpace().entitySet();
277 std::size_t intersectionIndex = 0;
278 for (const auto& is : intersections(lfsu.gridFunctionSpace().gridView(), eg.entity()))
279 {
280 Dune::PDELab::IntersectionGeometry<std::decay_t<decltype(is)>> ig(is, intersectionIndex++);
281 auto intersectionData = classifyIntersection(entitySet, is);
282 auto intersectionType = std::get<0>(intersectionData);
283
284 // Assemble the intersection part
285 switch (intersectionType){
287 if (LocalOperator::isLinear)
288 Dune::PDELab::impl::jacobianApplySkeleton(localOperator, ig, lfsu, z, lfsv, lfsu, z, lfsv, y, y);
289 else
290 Dune::PDELab::impl::jacobianApplySkeleton(localOperator, ig, lfsu, x, z, lfsv, lfsu, x, z, lfsv, y, y);
291 break;
293 if (LocalOperator::isLinear)
294 Dune::PDELab::impl::jacobianApplyBoundary(localOperator, ig, lfsu, z, lfsv, y);
295 else
296 Dune::PDELab::impl::jacobianApplyBoundary(localOperator, ig, lfsu, x, z, lfsv, y);
297 break;
298 default:
299 break;
300 }
301 }
302 }
303
304 } // namespace PDELab
305} // namespace Dune
306
307#endif
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
@ boundary
domain boundary intersection (neighbor() == false && boundary() == true)
@ skeleton
skeleton intersection (neighbor() == true && boundary() == false)
void assembleLocalDiagonalBlock(const LocalOperator &localOperator, const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, MAT &mat)
A function for assembling a single diagonal block.
Definition: blockdiagonalwrapper.hh:226
void applyLocalDiagonalBlock(const LocalOperator &localOperator, const EG &eg, const LFSU &lfsu, const X &x, const Z &z, const LFSV &lfsv, Y &y)
A function for applying a single diagonal block.
Definition: blockdiagonalwrapper.hh:261
std::tuple< IntersectionType, typename EntitySet::Element > classifyIntersection(const EntitySet &entity_set, const Intersection &is)
Classifies the type of an intersection wrt to the passed EntitySet.
Definition: intersectiontype.hh:37
std::enable_if_t< LOP::isLinear > jacobianApplyVolume(const LOP &lop, const EG &eg, const LFSU &lfsu, const X &z, const LFSV &lfsv, Y &y)
Definition: jacobianapplyhelper.hh:23
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
std::enable_if_t< LOP::isLinear > jacobianApplyBoundary(const LOP &lop, const IG &ig, const LFSU &lfsu_s, const X &z_s, const LFSV &lfsv_s, Y &y_s)
Definition: jacobianapplyhelper.hh:109
Wrap intersection.
Definition: geometrywrapper.hh:57
Definition: blockdiagonalwrapper.hh:19
auto data()
Definition: blockdiagonalwrapper.hh:31
const auto data() const
Definition: blockdiagonalwrapper.hh:36
const auto weight()
Definition: blockdiagonalwrapper.hh:26
void accumulate(const LFSU &lfsu, I i, const LFSV &lfsv, J j, Value value)
Definition: blockdiagonalwrapper.hh:50
BlockDiagonalAccumulationViewWrapper(AccumulationView &view, bool diagonal)
Definition: blockdiagonalwrapper.hh:22
void accumulate(const LFSU &lfsu, I i, Value value)
Definition: blockdiagonalwrapper.hh:42
A local operator that accumulates the block diagonal.
Definition: blockdiagonalwrapper.hh:82
void jacobian_apply_boundary(const IG &ig, const LFSU &lfsu_s, const X &z_s, const LFSV &lfsv_s, Y &y_s) const
Definition: blockdiagonalwrapper.hh:202
static constexpr bool doAlphaBoundary
Definition: blockdiagonalwrapper.hh:92
BlockDiagonalLocalOperatorWrapper(const LocalOperator &localOperator)
Construct new instance of class.
Definition: blockdiagonalwrapper.hh:108
BlockDiagonalLocalOperatorWrapper(const BlockDiagonalLocalOperatorWrapper &other)
Copy constructor.
Definition: blockdiagonalwrapper.hh:113
static constexpr bool doSkeletonTwoSided
Definition: blockdiagonalwrapper.hh:102
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, MAT &mat_ss) const
Definition: blockdiagonalwrapper.hh:152
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: blockdiagonalwrapper.hh:136
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const Z &z, const LFSV &lfsv, Y &y) const
Definition: blockdiagonalwrapper.hh:166
static constexpr bool doAlphaVolume
Definition: blockdiagonalwrapper.hh:90
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, MAT &mat) const
Definition: blockdiagonalwrapper.hh:126
static constexpr bool isLinear
Definition: blockdiagonalwrapper.hh:95
static constexpr bool doPatternVolume
Definition: blockdiagonalwrapper.hh:86
void jacobian_apply_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const Z &z_s, const LFSV &lfsv_s, Y &y_s) const
Definition: blockdiagonalwrapper.hh:210
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: blockdiagonalwrapper.hh:173
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: blockdiagonalwrapper.hh:187
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &z, const LFSV &lfsv, Y &y) const
Definition: blockdiagonalwrapper.hh:161
static constexpr bool doAlphaSkeleton
Definition: blockdiagonalwrapper.hh:91
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: blockdiagonalwrapper.hh:119
Default flags for all local operators.
Definition: flags.hh:19
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139