dune-pdelab  2.7-git
pointdiagonalwrapper.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_POINTDIAGONALWRAPPER_HH
4 #define DUNE_PDELAB_LOCALOPERATOR_POINTDIAGONALWRAPPER_HH
5 
6 namespace Dune {
7  namespace PDELab {
8 
9  namespace impl {
10  // This wraps a matrix accumulation view and only accumulates if diagonal
11  // is set to true and if the indices i and j are equal. Can be used to
12  // accumulate the point diagonal.
13  template <typename AccumulationView>
15  {
16  public:
17  PointDiagonalAccumulationViewWrapper(AccumulationView& view, bool diagonal)
18  : _view(view), _diagonal(diagonal)
19  {}
20 
21  template <typename LFSU, typename I, typename LFSV, typename J, typename Value>
22  void accumulate(const LFSU& lfsu, I i, const LFSV& lfsv, J j, Value value)
23  {
24  if (_diagonal && i == j){
25  _view.accumulate(lfsu, i, value);
26  }
27  }
28 
29  private:
30  AccumulationView& _view;
31  bool _diagonal;
32  };
33 
34  } // namespace impl
35 
36 
51  template <class LocalOperator>
54  {
55  public:
56  // We only want to assemble the point diagonal so we only need volume
57  // pattern.
58  static constexpr bool doPatternVolume = true;
59 
60  // For assembling the point diagonal correctly we might need to call
61  // volume, skeleton and boundary
62  static constexpr bool doAlphaVolume = LocalOperator::doAlphaVolume;
63  static constexpr bool doAlphaSkeleton = LocalOperator::doAlphaSkeleton;
64  static constexpr bool doAlphaBoundary = LocalOperator::doAlphaBoundary;
65 
66  // If the underlying lop is linear, this one will be linear too
67  static constexpr bool isLinear = LocalOperator::isLinear;
68 
69  // This wrapper is designed to use two sided assembly. If it was just
70  // about assembling the whole diagonal block matrix one sided assembly
71  // would be more efficient. For the implementation of matrix-free
72  // preconditioner we need to assemble a diagonal block locally for a
73  // given cell so we need to assemble in a two sided fashion.
74  static constexpr bool doSkeletonTwoSided = true;
75 
80  PointDiagonalLocalOperatorWrapper(const LocalOperator& localOperator)
81  : _localOperator(localOperator)
82 
83  {}
84 
87  : _localOperator(other._localOperator)
88  {}
89 
90  // define sparsity pattern of operator representation
91  template<typename LFSU, typename LFSV, typename LocalPattern>
92  void pattern_volume (const LFSU& lfsu, const LFSV& lfsv,
93  LocalPattern& pattern) const
94  {
95  _localOperator.pattern_volume(lfsu, lfsv, pattern);
96  }
97 
98  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
99  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
100  {
102  _localOperator.jacobian_volume(eg, lfsu, x, lfsv, view);
103  }
104 
105  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
106  void alpha_skeleton (const IG& ig,
107  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
108  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
109  R& r_s, R& r_n) const
110  {
112  impl::PointDiagonalAccumulationViewWrapper<R> view_other(r_s, false);
113  _localOperator.jacobian_skeleton(ig,
114  lfsu_s, x_s, lfsv_s,
115  lfsu_n, x_n, lfsv_n,
116  view_ss, view_other, view_other, view_other);
117  }
118 
119  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
120  void alpha_boundary (const IG& ig,
121  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
122  R& r_s) const
123  {
125  _localOperator.jacobian_boundary(ig, lfsu_s, x_s, lfsv_s, view);
126  }
127 
128  private:
130  const LocalOperator& _localOperator;
131  };
132 
138  template <typename LocalOperator, typename EG, typename LFSU, typename X, typename LFSV, typename Y>
139  void assembleLocalPointDiagonal(const LocalOperator& localOperator,
140  const EG& eg,
141  const LFSU& lfsu,
142  const X& x,
143  const LFSV& lfsv,
144  Y& y)
145  {
146  // Assemble the volume part
147  localOperator.alpha_volume(eg, lfsu, x, lfsv, y);
148 
149  // Iterate over the intersections
150  auto entitySet = lfsu.gridFunctionSpace().entitySet();
151  std::size_t intersectionIndex = 0;
152  for (const auto& is : intersections(lfsu.gridFunctionSpace().gridView(), eg.entity()))
153  {
155  typename std::remove_reference<
156  typename std::remove_const<
157  decltype(is)
158  >::type
159  >::type
160  > ig(is, intersectionIndex++);
161  auto intersectionData = classifyIntersection(entitySet, is);
162  auto intersectionType = std::get<0>(intersectionData);
163 
164  // Assemble the intersection part
165  switch (intersectionType){
167  localOperator.alpha_skeleton(ig, lfsu, x, lfsv, lfsu, x, lfsv, y, y);
168  break;
170  localOperator.alpha_boundary(ig, lfsu, x, lfsv, y);
171  break;
172  default:
173  break;
174  }
175  }
176  }
177 
178  } // namespace PDELab
179 } // namespace Dune
180 
181 #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)
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 assembleLocalPointDiagonal(const LocalOperator &localOperator, const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y)
A function for assembling the point diagonal of a single block.
Definition: pointdiagonalwrapper.hh:139
Wrap intersection.
Definition: geometrywrapper.hh:57
Default flags for all local operators.
Definition: flags.hh:19
Definition: pointdiagonalwrapper.hh:15
void accumulate(const LFSU &lfsu, I i, const LFSV &lfsv, J j, Value value)
Definition: pointdiagonalwrapper.hh:22
PointDiagonalAccumulationViewWrapper(AccumulationView &view, bool diagonal)
Definition: pointdiagonalwrapper.hh:17
A local operator that accumulates the point diagonal.
Definition: pointdiagonalwrapper.hh:54
void alpha_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, R &r_s, R &r_n) const
Definition: pointdiagonalwrapper.hh:106
PointDiagonalLocalOperatorWrapper(const LocalOperator &localOperator)
Construct new instance of class.
Definition: pointdiagonalwrapper.hh:80
static constexpr bool isLinear
Definition: pointdiagonalwrapper.hh:67
static constexpr bool doAlphaVolume
Definition: pointdiagonalwrapper.hh:62
static constexpr bool doPatternVolume
Definition: pointdiagonalwrapper.hh:58
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: pointdiagonalwrapper.hh:99
PointDiagonalLocalOperatorWrapper(const PointDiagonalLocalOperatorWrapper &other)
Copy constructor.
Definition: pointdiagonalwrapper.hh:86
static constexpr bool doAlphaBoundary
Definition: pointdiagonalwrapper.hh:64
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: pointdiagonalwrapper.hh:92
static constexpr bool doAlphaSkeleton
Definition: pointdiagonalwrapper.hh:63
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: pointdiagonalwrapper.hh:120
static constexpr bool doSkeletonTwoSided
Definition: pointdiagonalwrapper.hh:74
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139