diff --git a/dune/fufem/estimators/hierarchicestimatorbase.hh b/dune/fufem/estimators/hierarchicestimatorbase.hh index bfae7c5877037c7100d9c64863226b19c52bf48b..943cabf84ac7c875cde9413ce85171aba47713b8 100644 --- a/dune/fufem/estimators/hierarchicestimatorbase.hh +++ b/dune/fufem/estimators/hierarchicestimatorbase.hh @@ -12,38 +12,31 @@ #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh> #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh> - template <class EnlargedGlobalBasis, int block_size, class field_type> class HierarchicEstimatorBase { - typedef typename EnlargedGlobalBasis::GridView GridView; + using GridView = typename EnlargedGlobalBasis::GridView; + using Grid = typename GridView::Grid; + using Lfe = typename EnlargedGlobalBasis::LocalFiniteElement; - typedef typename GridView::template Codim<0>::Iterator ElementIterator; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::IntersectionIterator IntersectionIterator; enum {dim = GridView::dimension}; - enum {dimworld = GridView::dimensionworld}; - typedef typename GridView::ctype CType; - - typedef Dune::FieldMatrix<field_type,block_size,block_size> BlockType; - typedef Dune::FieldVector<field_type,block_size> FieldVectorType; + using CType = typename GridView::ctype; + using BlockType = Dune::FieldMatrix<field_type,block_size,block_size>; + using FieldVectorType = Dune::FieldVector<field_type,block_size>; public: static void preconditionedDefectProblem(const EnlargedGlobalBasis& enlargedGlobalBasis, - const Dune::BlockVector<Dune::FieldVector<field_type,block_size> >& x, + const Dune::BlockVector<FieldVectorType>& x, const Dune::VirtualFunction<Dune::FieldVector<CType,dimworld>, Dune::FieldVector<field_type,block_size> >* volumeTerm, const Dune::VirtualFunction<Dune::FieldVector<CType,dimworld>, Dune::FieldVector<field_type,block_size> >* neumannTerm, - Dune::BlockVector<Dune::FieldVector<field_type,block_size> >& residual, + Dune::BlockVector<FieldVectorType>& residual, Dune::BDMatrix<BlockType>& matrix, - LocalOperatorAssembler<typename GridView::Grid, - typename EnlargedGlobalBasis::LocalFiniteElement, - typename EnlargedGlobalBasis::LocalFiniteElement, - Dune::FieldMatrix<field_type,block_size,block_size> >* localStiffness) + LocalOperatorAssembler<Grid, Lfe, Lfe, BlockType>* localStiffness) { const GridView& gridView = enlargedGlobalBasis.getGridView(); @@ -56,20 +49,17 @@ public: matrix = 0; // local assembler for the volume term - L2FunctionalAssembler<typename GridView::Grid, typename EnlargedGlobalBasis::LocalFiniteElement, Dune::FieldVector<field_type,block_size> > localL2FunctionalAssembler(*volumeTerm); + L2FunctionalAssembler<Grid, Lfe, FieldVectorType> localL2FunctionalAssembler(*volumeTerm); // //////////////////////////////////////////////////////////////////////// // Loop over all elements and assemble diagonal entries and residual // //////////////////////////////////////////////////////////////////////// - ElementIterator eIt = gridView.template begin<0>(); - ElementIterator eEndIt = gridView.template end<0>(); - - for (; eIt!=eEndIt; ++eIt) { + for (const auto& e : elements(gridView)) { // Get set of shape functions - const typename EnlargedGlobalBasis::LocalFiniteElement& localFiniteElement = enlargedGlobalBasis.getLocalFiniteElement(*eIt); - const typename EnlargedGlobalBasis::LocalFiniteElement::Traits::LocalBasisType& localBasis = localFiniteElement.localBasis(); + const auto& localFiniteElement = enlargedGlobalBasis.getLocalFiniteElement(e); + const auto& localBasis = localFiniteElement.localBasis(); // ///////////////////////////////////////////////////////////////////////// // Assemble right hand side in the enlarged space @@ -79,14 +69,14 @@ public: // vector for the local rhs from the volume term //Dune::BlockVector<Dune::FieldVector<field_type,block_size> > localVolumeVector(localBasis.size()); - typename L2FunctionalAssembler<typename GridView::Grid, EnlargedGlobalBasis, FieldVectorType>::LocalVector localVolumeVector(localBasis.size()); + typename L2FunctionalAssembler<Grid, EnlargedGlobalBasis, FieldVectorType>::LocalVector localVolumeVector(localBasis.size()); // assemble local rhs from the volume term - localL2FunctionalAssembler.assemble(*eIt, localVolumeVector, localFiniteElement); + localL2FunctionalAssembler.assemble(e, localVolumeVector, localFiniteElement); // add to global vector for (size_t i=0; i<localBasis.size(); i++) { - int globalRow = enlargedGlobalBasis.index(*eIt, i); + int globalRow = enlargedGlobalBasis.index(e, i); residual[globalRow] += localVolumeVector[i]; } @@ -95,20 +85,20 @@ public: if (neumannTerm) { const int quadOrder = 6; - NeumannBoundaryAssembler<typename GridView::Grid,FieldVectorType> neumannAssembler(*neumannTerm,quadOrder); + NeumannBoundaryAssembler<Grid, FieldVectorType> neumannAssembler(*neumannTerm,quadOrder); - auto nIt = gridView.ibegin(*eIt); - auto nEndIt = gridView.iend(*eIt); + auto nIt = gridView.ibegin(e); + auto nEndIt = gridView.iend(e); for (;nIt != nEndIt; ++nIt) { if (!nIt->boundary()) continue; - typename NeumannBoundaryAssembler<typename GridView::Grid, FieldVectorType>::LocalVector localNeumannVector(localBasis.size()); + typename NeumannBoundaryAssembler<Grid, FieldVectorType>::LocalVector localNeumannVector(localBasis.size()); neumannAssembler.assemble(nIt,localNeumannVector,localFiniteElement); for (size_t i=0; i<localBasis.size(); i++) { - int globalRow = enlargedGlobalBasis.index(*eIt, i); + int globalRow = enlargedGlobalBasis.index(e, i); residual[globalRow] += localNeumannVector[i]; } } @@ -123,24 +113,23 @@ public: Dune::BlockVector<FieldVectorType> localSolution(localBasis.size()); for (size_t i=0; i<localBasis.size(); i++) - localSolution[i] = x[enlargedGlobalBasis.index(*eIt, i)]; + localSolution[i] = x[enlargedGlobalBasis.index(e, i)]; // Assemble element matrix with quadratic Lagrangian elements Dune::Matrix<BlockType> localMatrix(localBasis.size(),localBasis.size()); - localStiffness->assemble(*eIt,localSolution, localMatrix, + localStiffness->assemble(e,localSolution, localMatrix, localFiniteElement, localFiniteElement); // Add to residual vector for (size_t i=0; i<localBasis.size(); i++) { - int globalRow = enlargedGlobalBasis.index(*eIt, i); + int globalRow = enlargedGlobalBasis.index(e, i); for (size_t j=0; j<localBasis.size(); j++) { - int globalCol = enlargedGlobalBasis.index(*eIt, j); + int globalCol = enlargedGlobalBasis.index(e, j); - //residual[globalRow] -= localStiffness.mat(i,j) * x[globalCol]; localMatrix[i][j].mmv(x[globalCol], residual[globalRow]); } @@ -175,13 +164,10 @@ public: // Loop over all elements and assemble diagonal entries // //////////////////////////////////////////////////////////////////////// - ElementIterator eIt = gridView.template begin<0>(); - ElementIterator eEndIt = gridView.template end<0>(); - - for (; eIt!=eEndIt; ++eIt) { + for (const auto& e : elements(gridView)) { // Get set of shape functions - const typename EnlargedGlobalBasis::LocalFiniteElement& enlargedGlobalFiniteElement = enlargedGlobalBasis.getLocalFiniteElement(*eIt); + const auto& enlargedGlobalFiniteElement = enlargedGlobalBasis.getLocalFiniteElement(e); localMatrix.setSize(enlargedGlobalFiniteElement.localBasis().size(), enlargedGlobalFiniteElement.localBasis().size()); @@ -189,17 +175,17 @@ public: Dune::BlockVector<FieldVectorType> localSolution(enlargedGlobalFiniteElement.localBasis().size()); for (size_t i=0; i<enlargedGlobalFiniteElement.localBasis().size(); i++) - localSolution[i] = x[enlargedGlobalBasis.index(*eIt,i)]; + localSolution[i] = x[enlargedGlobalBasis.index(e, i)]; // Assemble element matrix with quadratic Lagrangian elements - localStiffness.assemble(*eIt, localSolution, localMatrix, + localStiffness.assemble(e, localSolution, localMatrix, enlargedGlobalFiniteElement, enlargedGlobalFiniteElement); // Add to the total energy for (size_t i=0; i<enlargedGlobalFiniteElement.localBasis().size(); i++) { - int globalRow = enlargedGlobalBasis.index(*eIt, i); + int globalRow = enlargedGlobalBasis.index(e, i); FieldVectorType smallTmp(0); localMatrix[i][i].umv(x[globalRow], smallTmp); energyNormSquared += smallTmp * x[globalRow];