dune-pdelab  2.7-git
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 
8 namespace 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:
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.
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
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
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
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