4#ifndef DUNE_PDELAB_LOCALOPERATOR_BLOCKDIAGONALWRAPPER_HH
5#define DUNE_PDELAB_LOCALOPERATOR_BLOCKDIAGONALWRAPPER_HH
7#include <dune/pdelab/common/intersectiontype.hh>
8#include <dune/pdelab/localoperator/jacobianapplyhelper.hh>
17 template <
typename AccumulationView>
18 class BlockDiagonalAccumulationViewWrapper
22 BlockDiagonalAccumulationViewWrapper(AccumulationView& view,
bool diagonal)
23 : _view(view), _diagonal(diagonal)
28 return _view.weight();
36 const auto data()
const
41 template <
typename LFSU,
typename I,
typename Value>
42 void accumulate(
const LFSU& lfsu, I i, Value value)
45 _view.accumulate(lfsu, i, value);
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)
53 _view.accumulate(lfsu, i, lfsv, j, value);
58 AccumulationView& _view;
79 template <
class LocalOperator>
86 static constexpr bool doPatternVolume =
true;
90 static constexpr bool doAlphaVolume = LocalOperator::doAlphaVolume;
91 static constexpr bool doAlphaSkeleton = LocalOperator::doAlphaSkeleton;
92 static constexpr bool doAlphaBoundary = LocalOperator::doAlphaBoundary;
95 static constexpr bool isLinear = LocalOperator::isLinear;
102 static constexpr bool doSkeletonTwoSided =
true;
109 : _localOperator(localOperator)
114 : _localOperator(other._localOperator)
118 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
119 void pattern_volume (
const LFSU& lfsu,
const LFSV& lfsv,
120 LocalPattern& pattern)
const
122 _localOperator.pattern_volume(lfsu, lfsv, pattern);
125 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename MAT>
126 void jacobian_volume (
const EG& eg,
132 _localOperator.jacobian_volume(eg,lfsu,x,lfsv,mat);
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
146 impl::BlockDiagonalAccumulationViewWrapper<MAT> view_ss(mat_ss,
true);
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);
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,
156 _localOperator.jacobian_boundary(ig, lfsu_s, x_s, lfsv_s, mat_ss);
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
163 Dune::PDELab::impl::jacobianApplyVolume(_localOperator, eg, lfsu, z, lfsv, y);
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
168 Dune::PDELab::impl::jacobianApplyVolume(_localOperator, eg, lfsu, x, z, lfsv, y);
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
178 impl::BlockDiagonalAccumulationViewWrapper<Y> view_s(y_s,
true);
179 impl::BlockDiagonalAccumulationViewWrapper<Y> view_n(y_n,
false);
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);
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
192 impl::BlockDiagonalAccumulationViewWrapper<Y> view_s(y_s,
true);
193 impl::BlockDiagonalAccumulationViewWrapper<Y> view_n(y_n,
false);
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);
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,
206 Dune::PDELab::impl::jacobianApplyBoundary(_localOperator, ig, lfsu_s, z_s, lfsv_s, y_s);
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,
214 Dune::PDELab::impl::jacobianApplyBoundary(_localOperator, ig, lfsu_s, x_s, z_s, lfsv_s, y_s);
220 const LocalOperator& _localOperator;
225 template <
typename LocalOperator,
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename MAT>
226 void assembleLocalDiagonalBlock(
const LocalOperator& localOperator,
234 localOperator.jacobian_volume(eg, lfsu, x, lfsv, mat);
237 auto entitySet = lfsu.gridFunctionSpace().entitySet();
238 std::size_t intersectionIndex = 0;
239 for (
const auto& is : intersections(lfsu.gridFunctionSpace().gridView(), eg.entity()))
242 auto intersectionData = classifyIntersection(entitySet, is);
243 auto intersectionType = std::get<0>(intersectionData);
246 switch (intersectionType){
247 case Dune::PDELab::IntersectionType::skeleton:
248 localOperator.jacobian_skeleton(ig, lfsu, x, lfsv, lfsu, x, lfsv, mat, mat, mat, mat);
250 case Dune::PDELab::IntersectionType::boundary:
251 localOperator.jacobian_boundary(ig, lfsu, x, lfsv, mat);
260 template<
typename LocalOperator,
typename EG,
typename LFSU,
typename X,
typename Z,
typename LFSV,
typename Y>
261 void applyLocalDiagonalBlock(
const LocalOperator& localOperator,
270 if (LocalOperator::isLinear)
271 Dune::PDELab::impl::jacobianApplyVolume(localOperator, eg, lfsu, z, lfsv, y);
273 Dune::PDELab::impl::jacobianApplyVolume(localOperator, eg, lfsu, x, z, lfsv, y);
276 auto entitySet = lfsu.gridFunctionSpace().entitySet();
277 std::size_t intersectionIndex = 0;
278 for (
const auto& is : intersections(lfsu.gridFunctionSpace().gridView(), eg.entity()))
281 auto intersectionData = classifyIntersection(entitySet, is);
282 auto intersectionType = std::get<0>(intersectionData);
285 switch (intersectionType){
286 case Dune::PDELab::IntersectionType::skeleton:
287 if (LocalOperator::isLinear)
288 Dune::PDELab::impl::jacobianApplySkeleton(localOperator, ig, lfsu, z, lfsv, lfsu, z, lfsv, y, y);
290 Dune::PDELab::impl::jacobianApplySkeleton(localOperator, ig, lfsu, x, z, lfsv, lfsu, x, z, lfsv, y, y);
292 case Dune::PDELab::IntersectionType::boundary:
293 if (LocalOperator::isLinear)
294 Dune::PDELab::impl::jacobianApplyBoundary(localOperator, ig, lfsu, z, lfsv, y);
296 Dune::PDELab::impl::jacobianApplyBoundary(localOperator, ig, lfsu, x, z, lfsv, y);
A local operator that accumulates the block diagonal.
Definition: blockdiagonalwrapper.hh:82
BlockDiagonalLocalOperatorWrapper(const LocalOperator &localOperator)
Construct new instance of class.
Definition: blockdiagonalwrapper.hh:108
BlockDiagonalLocalOperatorWrapper(const BlockDiagonalLocalOperatorWrapper &other)
Copy constructor.
Definition: blockdiagonalwrapper.hh:113
Wrap intersection.
Definition: geometrywrapper.hh:57
Default flags for all local operators.
Definition: flags.hh:19
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
Dune namespace.
Definition: alignedallocator.hh:11