dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
6namespace 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 {
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