Skip to content
Snippets Groups Projects
Commit f63cbffb authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Cleanup a bit using range-based for-loops

parent b6ed3dc9
No related branches found
No related tags found
No related merge requests found
...@@ -12,38 +12,31 @@ ...@@ -12,38 +12,31 @@
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh> #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh> #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
template <class EnlargedGlobalBasis, int block_size, class field_type> template <class EnlargedGlobalBasis, int block_size, class field_type>
class HierarchicEstimatorBase 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 {dim = GridView::dimension};
enum {dimworld = GridView::dimensionworld}; enum {dimworld = GridView::dimensionworld};
typedef typename GridView::ctype CType; using CType = typename GridView::ctype;
using BlockType = Dune::FieldMatrix<field_type,block_size,block_size>;
typedef Dune::FieldMatrix<field_type,block_size,block_size> BlockType; using FieldVectorType = Dune::FieldVector<field_type,block_size>;
typedef Dune::FieldVector<field_type,block_size> FieldVectorType;
public: public:
static void preconditionedDefectProblem(const EnlargedGlobalBasis& enlargedGlobalBasis, 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> >* volumeTerm,
const Dune::VirtualFunction<Dune::FieldVector<CType,dimworld>, Dune::FieldVector<field_type,block_size> >* neumannTerm, 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, Dune::BDMatrix<BlockType>& matrix,
LocalOperatorAssembler<typename GridView::Grid, LocalOperatorAssembler<Grid, Lfe, Lfe, BlockType>* localStiffness)
typename EnlargedGlobalBasis::LocalFiniteElement,
typename EnlargedGlobalBasis::LocalFiniteElement,
Dune::FieldMatrix<field_type,block_size,block_size> >* localStiffness)
{ {
const GridView& gridView = enlargedGlobalBasis.getGridView(); const GridView& gridView = enlargedGlobalBasis.getGridView();
...@@ -56,20 +49,17 @@ public: ...@@ -56,20 +49,17 @@ public:
matrix = 0; matrix = 0;
// local assembler for the volume term // 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 // Loop over all elements and assemble diagonal entries and residual
// //////////////////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////////////////
ElementIterator eIt = gridView.template begin<0>(); for (const auto& e : elements(gridView)) {
ElementIterator eEndIt = gridView.template end<0>();
for (; eIt!=eEndIt; ++eIt) {
// Get set of shape functions // Get set of shape functions
const typename EnlargedGlobalBasis::LocalFiniteElement& localFiniteElement = enlargedGlobalBasis.getLocalFiniteElement(*eIt); const auto& localFiniteElement = enlargedGlobalBasis.getLocalFiniteElement(e);
const typename EnlargedGlobalBasis::LocalFiniteElement::Traits::LocalBasisType& localBasis = localFiniteElement.localBasis(); const auto& localBasis = localFiniteElement.localBasis();
// ///////////////////////////////////////////////////////////////////////// // /////////////////////////////////////////////////////////////////////////
// Assemble right hand side in the enlarged space // Assemble right hand side in the enlarged space
...@@ -79,14 +69,14 @@ public: ...@@ -79,14 +69,14 @@ public:
// vector for the local rhs from the volume term // vector for the local rhs from the volume term
//Dune::BlockVector<Dune::FieldVector<field_type,block_size> > localVolumeVector(localBasis.size()); //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 // assemble local rhs from the volume term
localL2FunctionalAssembler.assemble(*eIt, localVolumeVector, localFiniteElement); localL2FunctionalAssembler.assemble(e, localVolumeVector, localFiniteElement);
// add to global vector // add to global vector
for (size_t i=0; i<localBasis.size(); i++) { 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]; residual[globalRow] += localVolumeVector[i];
} }
...@@ -95,20 +85,20 @@ public: ...@@ -95,20 +85,20 @@ public:
if (neumannTerm) { if (neumannTerm) {
const int quadOrder = 6; const int quadOrder = 6;
NeumannBoundaryAssembler<typename GridView::Grid,FieldVectorType> neumannAssembler(*neumannTerm,quadOrder); NeumannBoundaryAssembler<Grid, FieldVectorType> neumannAssembler(*neumannTerm,quadOrder);
auto nIt = gridView.ibegin(*eIt); auto nIt = gridView.ibegin(e);
auto nEndIt = gridView.iend(*eIt); auto nEndIt = gridView.iend(e);
for (;nIt != nEndIt; ++nIt) { for (;nIt != nEndIt; ++nIt) {
if (!nIt->boundary()) if (!nIt->boundary())
continue; continue;
typename NeumannBoundaryAssembler<typename GridView::Grid, FieldVectorType>::LocalVector localNeumannVector(localBasis.size()); typename NeumannBoundaryAssembler<Grid, FieldVectorType>::LocalVector localNeumannVector(localBasis.size());
neumannAssembler.assemble(nIt,localNeumannVector,localFiniteElement); neumannAssembler.assemble(nIt,localNeumannVector,localFiniteElement);
for (size_t i=0; i<localBasis.size(); i++) { 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]; residual[globalRow] += localNeumannVector[i];
} }
} }
...@@ -123,24 +113,23 @@ public: ...@@ -123,24 +113,23 @@ public:
Dune::BlockVector<FieldVectorType> localSolution(localBasis.size()); Dune::BlockVector<FieldVectorType> localSolution(localBasis.size());
for (size_t i=0; i<localBasis.size(); i++) 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 // Assemble element matrix with quadratic Lagrangian elements
Dune::Matrix<BlockType> localMatrix(localBasis.size(),localBasis.size()); Dune::Matrix<BlockType> localMatrix(localBasis.size(),localBasis.size());
localStiffness->assemble(*eIt,localSolution, localMatrix, localStiffness->assemble(e,localSolution, localMatrix,
localFiniteElement, localFiniteElement,
localFiniteElement); localFiniteElement);
// Add to residual vector // Add to residual vector
for (size_t i=0; i<localBasis.size(); i++) { 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++) { 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]); localMatrix[i][j].mmv(x[globalCol], residual[globalRow]);
} }
...@@ -175,13 +164,10 @@ public: ...@@ -175,13 +164,10 @@ public:
// Loop over all elements and assemble diagonal entries // Loop over all elements and assemble diagonal entries
// //////////////////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////////////////
ElementIterator eIt = gridView.template begin<0>(); for (const auto& e : elements(gridView)) {
ElementIterator eEndIt = gridView.template end<0>();
for (; eIt!=eEndIt; ++eIt) {
// Get set of shape functions // 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()); localMatrix.setSize(enlargedGlobalFiniteElement.localBasis().size(), enlargedGlobalFiniteElement.localBasis().size());
...@@ -189,17 +175,17 @@ public: ...@@ -189,17 +175,17 @@ public:
Dune::BlockVector<FieldVectorType> localSolution(enlargedGlobalFiniteElement.localBasis().size()); Dune::BlockVector<FieldVectorType> localSolution(enlargedGlobalFiniteElement.localBasis().size());
for (size_t i=0; i<enlargedGlobalFiniteElement.localBasis().size(); i++) 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 // Assemble element matrix with quadratic Lagrangian elements
localStiffness.assemble(*eIt, localSolution, localMatrix, localStiffness.assemble(e, localSolution, localMatrix,
enlargedGlobalFiniteElement, enlargedGlobalFiniteElement,
enlargedGlobalFiniteElement); enlargedGlobalFiniteElement);
// Add to the total energy // Add to the total energy
for (size_t i=0; i<enlargedGlobalFiniteElement.localBasis().size(); i++) { 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); FieldVectorType smallTmp(0);
localMatrix[i][i].umv(x[globalRow], smallTmp); localMatrix[i][i].umv(x[globalRow], smallTmp);
energyNormSquared += smallTmp * x[globalRow]; energyNormSquared += smallTmp * x[globalRow];
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment