dune-pdelab  2.7-git
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 
10 namespace 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>
173  void jacobian_apply_skeleton(const IG& ig,
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>
187  void jacobian_apply_skeleton(const IG& ig,
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>
202  void jacobian_apply_boundary (const IG& ig,
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>
210  void jacobian_apply_boundary (const IG& ig,
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
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
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::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 > 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 > 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