3#ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ASSEMBLEDBLOCKJACOBIPREONDITIONER_HH
4#define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ASSEMBLEDBLOCKJACOBIPREONDITIONER_HH
10#include <dune/pdelab/localoperator/blockdiagonalwrapper.hh>
22 class WeightedEigenMatrixAccumulationView
28 using size_type =
typename Container::StorageIndex;
29 using weight_type = value_type;
31 const weight_type& weight()
38 return _container.data();
43 return _container.data();
51 template <
typename LFSV,
typename LFSU>
52 void accumulate(
const LFSV& lfsv, size_type i,
const LFSU& lfsu, size_type j, value_type value)
54 _container(i,j) += _weight * value;
57 WeightedEigenMatrixAccumulationView(Container& container, weight_type weight)
58 : _container(container)
63 Container& _container;
96 template<
class BlockDiagonalLocalOperator,
class Gr
idFunctionSpace,
class DomainField>
97 class AssembledBlockJacobiPreconditionerLocalOperator
102 using GridView =
typename GridFunctionSpace::Traits::GridViewType;
103 using Grid =
typename GridView::Traits::Grid;
104 using EntityType =
typename Grid::template Codim<0>::Entity;
105 using value_type = DomainField;
106 using DiagonalBlock = ::Eigen::Matrix<value_type, ::Eigen::Dynamic, ::Eigen::Dynamic, ::Eigen::RowMajor>;
107 using PermutationIndices = ::Eigen::Matrix<int, ::Eigen::Dynamic, 1>;
108 using TupleType = std::tuple<DiagonalBlock, PermutationIndices>;
116 static constexpr bool doPatternVolume =
true;
117 static constexpr bool doAlphaVolume =
true;
118 static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
130 AssembledBlockJacobiPreconditionerLocalOperator(
const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
131 const GridFunctionSpace& gridFunctionSpace,
132 const bool verbose=0)
133 : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
134 , _gridFunctionSpace(gridFunctionSpace)
135 , _mapper(gridFunctionSpace.gridView())
136 , _precCache(_mapper.
size())
138 , _requireSetup(true)
149 return _requireSetup;
151 void setRequireSetup(
bool requireSetup)
153 _requireSetup = requireSetup;
162 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
163 void alpha_volume(
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, Y& y)
const
166 const auto size = lfsu.size();
167 assert(lfsv.size() ==
size);
170 DiagonalBlock diagonalBlock = DiagonalBlock::Constant(
size,
size, 0.0);
171 impl::WeightedEigenMatrixAccumulationView<DiagonalBlock> viewDiagonalBlock(diagonalBlock, y.weight());
172 assembleLocalDiagonalBlock(_blockDiagonalLocalOperator, eg, lfsu, x, lfsv, viewDiagonalBlock);
175 ::Eigen::PartialPivLU<::Eigen::Ref<DiagonalBlock>> lu(diagonalBlock);
176 auto inversePermutationIndices = lu.permutationP().indices();
179 PermutationIndices permutationIndices(
size);
180 for(
int i=0; i<inversePermutationIndices.size(); ++i)
181 permutationIndices[inversePermutationIndices[i]] = i;
184 std::size_t cache_idx = _mapper.index(eg.entity());
185 _precCache[cache_idx] = std::make_tuple(diagonalBlock, permutationIndices);
188 template<
typename EG,
typename LFSU,
typename X,
typename Z,
typename LFSV,
typename Y>
189 void jacobian_apply_volume(
const EG& eg,
const LFSU& lfsu,
const X& x,
const Z& z,
const LFSV& lfsv, Y& y)
const
193 jacobian_apply_volume(eg,lfsu,z,lfsv,y);
202 template<
typename EG,
typename LFSU,
typename Z,
typename LFSV,
typename Y>
203 void jacobian_apply_volume(
const EG& eg,
const LFSU& lfsu,
const Z& z,
const LFSV& lfsv, Y& y)
const
205 assert(not _requireSetup);
208 std::size_t cache_idx = _mapper.index(eg.entity());
209 const auto& cacheEntry = _precCache[cache_idx];
210 auto diagonalBlockLU = std::get<0>(cacheEntry);
211 auto permutationIndices = std::get<1>(cacheEntry);
214 const int size = lfsu.size();
215 const auto& z_container = accessBaseContainer(z);
216 std::vector<value_type> ztmp_container(
size);
217 auto& y_container = accessBaseContainer(y);
220 for(
int i=0; i<
size; ++i){
221 value_type rhs(z_container[permutationIndices[i]]);
222 for(
int j=0; j<i; ++j)
223 rhs -= diagonalBlockLU(i, j) * ztmp_container[j];
224 ztmp_container[i] = rhs;
227 for(
int i=
size-1; i>=0; i--){
228 value_type rhs(ztmp_container[i]);
229 for(
int j=
size-1; j>i; j--)
230 rhs -= diagonalBlockLU(i, j) * y_container[j];
231 y_container[i] = rhs/diagonalBlockLU(i,i);
236 BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
237 const GridFunctionSpace& _gridFunctionSpace;
239 mutable std::vector<TupleType> _precCache;
sparsity pattern generator
Definition: pattern.hh:14
Default flags for all local operators.
Definition: flags.hh:19
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Mapper classes are used to attach data to a grid.