Skip to content
Snippets Groups Projects
Commit f80d4a1b authored by Oliver Sander's avatar Oliver Sander
Browse files

[cleanup] Modernize the code with a bit of C++11

In particular:
- auto
- range-based for
parent 9b3530ea
No related branches found
No related tags found
No related merge requests found
...@@ -35,33 +35,27 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix) ...@@ -35,33 +35,27 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
{ {
static const int dim = GridView::Grid::dimension; static const int dim = GridView::Grid::dimension;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache; typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
typedef typename FiniteElementCache::FiniteElementType FiniteElement; typedef typename FiniteElementCache::FiniteElementType FiniteElement;
const IndexSet& indexSet = gridView.indexSet(); const auto& indexSet = gridView.indexSet();
FiniteElementCache cache; FiniteElementCache cache;
int size = indexSet.size(dim); int size = indexSet.size(dim);
Dune::MatrixIndexSet indices(size, size); Dune::MatrixIndexSet indices(size, size);
ElementIterator it = gridView.template begin<0>(); for (const auto& element : elements(gridView))
ElementIterator end = gridView.template end<0>();
for (; it != end; ++it)
{ {
const Element& e = *it; const FiniteElement& fe = cache.get(element.type());
const FiniteElement& fe = cache.get(e.type());
int localSize = fe.localBasis().size(); int localSize = fe.size();
for (int i = 0; i < localSize; ++i) for (int i = 0; i < localSize; ++i)
{ {
int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim); int iGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(i).subEntity(), dim);
for (int j = 0; j < localSize; ++j) for (int j = 0; j < localSize; ++j)
{ {
int jGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(j).subEntity(), dim); int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
indices.add(iGlobal, jGlobal); indices.add(iGlobal, jGlobal);
indices.add(jGlobal, iGlobal); indices.add(jGlobal, iGlobal);
} }
...@@ -77,39 +71,31 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix) ...@@ -77,39 +71,31 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
static const int dim = GridView::Grid::dimension; static const int dim = GridView::Grid::dimension;
static const int dimworld = GridView::Grid::dimensionworld; static const int dimworld = GridView::Grid::dimensionworld;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element; typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache; typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
typedef typename FiniteElementCache::FiniteElementType FiniteElement; typedef typename FiniteElementCache::FiniteElementType FiniteElement;
typedef typename Dune::QuadratureRule<double,dim> QuadratureRule; typedef typename Dune::FieldVector<double,dim> LocalCoordinate;
typedef typename Dune::template FieldVector<double,dim> LocalCoordinate; typedef typename Dune::FieldVector<double,dimworld> GlobalCoordinate;
typedef typename Dune::template FieldVector<double,dimworld> GlobalCoordinate;
typedef typename Element::Geometry::JacobianInverseTransposed JacobianInverseTransposed; typedef typename Element::Geometry::JacobianInverseTransposed JacobianInverseTransposed;
typedef typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType; typedef typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType;
const IndexSet& indexSet = gridView.indexSet(); const auto& indexSet = gridView.indexSet();
FiniteElementCache cache; FiniteElementCache cache;
ElementIterator it = gridView.template begin<0>(); for (const auto& element : elements(gridView))
ElementIterator end = gridView.template end<0>();
for (; it != end; ++it)
{ {
const Element& e = *it; const Geometry geometry = element.geometry();
const Geometry geometry = e.geometry(); const FiniteElement& fe = cache.get(element.type());
const FiniteElement& fe = cache.get(e.type());
int localSize = fe.localBasis().size(); int localSize = fe.size();
// get quadrature rule of appropiate order (P1/Q1) // get quadrature rule of appropiate order (P1/Q1)
int order = 0; int order = (element.type().isSimplex())
if (e.type().isSimplex()) ? 2*(1-1)
order = 2*(1-1); : 2*(dim-1);
else const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
order = 2*(dim-1);
const QuadratureRule& quad = Dune::template QuadratureRules<double, dim>::rule(e.type(), order);
// store gradients of shape functions and base functions // store gradients of shape functions and base functions
std::vector<JacobianType> referenceGradients(localSize); std::vector<JacobianType> referenceGradients(localSize);
...@@ -137,10 +123,10 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix) ...@@ -137,10 +123,10 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
double z = quad[pt].weight() * integrationElement; double z = quad[pt].weight() * integrationElement;
for (int i = 0; i < localSize; ++i) for (int i = 0; i < localSize; ++i)
{ {
int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim); int iGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(i).subEntity(), dim);
for (int j = i+1; j < localSize; ++j) for (int j = i+1; j < localSize; ++j)
{ {
int jGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(j).subEntity(), dim); int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
double zij = (gradients[i] * gradients[j]) * z; double zij = (gradients[i] * gradients[j]) * z;
matrix[iGlobal][jGlobal] += zij; matrix[iGlobal][jGlobal] += zij;
...@@ -159,37 +145,28 @@ void assemblePQ1Mass(const GridView& gridView, Matrix& matrix) ...@@ -159,37 +145,28 @@ void assemblePQ1Mass(const GridView& gridView, Matrix& matrix)
{ {
static const int dim = GridView::Grid::dimension; static const int dim = GridView::Grid::dimension;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache; typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
typedef typename FiniteElementCache::FiniteElementType FiniteElement; typedef typename FiniteElementCache::FiniteElementType FiniteElement;
typedef typename Dune::QuadratureRule<double,dim> QuadratureRule; typedef typename Dune::FieldVector<double,dim> LocalCoordinate;
typedef typename Dune::template FieldVector<double,dim> LocalCoordinate;
typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType; typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
const IndexSet& indexSet = gridView.indexSet(); const auto& indexSet = gridView.indexSet();
FiniteElementCache cache; FiniteElementCache cache;
ElementIterator it = gridView.template begin<0>(); for (const auto& element : elements(gridView))
ElementIterator end = gridView.template end<0>();
for (; it != end; ++it)
{ {
const Element& e = *it; const Geometry geometry = element.geometry();
const Geometry geometry = e.geometry(); const FiniteElement& fe = cache.get(element.type());
const FiniteElement& fe = cache.get(e.type());
int localSize = fe.localBasis().size(); int localSize = fe.size();
// get quadrature rule of appropiate order (P1/Q1) // get quadrature rule of appropriate order (P1/Q1)
int order = 0; int order = (element.type().isSimplex())
if (e.type().isSimplex()) ? 2*1
order = 2*1; : 2*dim;
else const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
order = 2*dim;
const QuadratureRule& quad = Dune::template QuadratureRules<double, dim>::rule(e.type(), order);
// store values of shape functions // store values of shape functions
std::vector<RangeType> values(localSize); std::vector<RangeType> values(localSize);
...@@ -209,11 +186,11 @@ void assemblePQ1Mass(const GridView& gridView, Matrix& matrix) ...@@ -209,11 +186,11 @@ void assemblePQ1Mass(const GridView& gridView, Matrix& matrix)
double z = quad[pt].weight() * integrationElement; double z = quad[pt].weight() * integrationElement;
for (int i = 0; i < localSize; ++i) for (int i = 0; i < localSize; ++i)
{ {
int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim); int iGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(i).subEntity(), dim);
double zi = values[i]*z; double zi = values[i]*z;
for (int j = i+1; j < localSize; ++j) for (int j = i+1; j < localSize; ++j)
{ {
int jGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(j).subEntity(), dim); int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
double zij = values[j] * zi; double zij = values[j] * zi;
matrix[iGlobal][jGlobal] += zij; matrix[iGlobal][jGlobal] += zij;
...@@ -232,38 +209,30 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f) ...@@ -232,38 +209,30 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
{ {
static const int dim = GridView::Grid::dimension; static const int dim = GridView::Grid::dimension;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache; typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
typedef typename FiniteElementCache::FiniteElementType FiniteElement; typedef typename FiniteElementCache::FiniteElementType FiniteElement;
typedef typename Dune::QuadratureRule<double,dim> QuadratureRule;
typedef typename Dune::template FieldVector<double,dim> LocalCoordinate; typedef typename Dune::template FieldVector<double,dim> LocalCoordinate;
typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType; typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
typedef typename Function::RangeType FunctionRangeType; typedef typename Function::RangeType FunctionRangeType;
const IndexSet& indexSet = gridView.indexSet(); const auto& indexSet = gridView.indexSet();
FiniteElementCache cache; FiniteElementCache cache;
ElementIterator it = gridView.template begin<0>(); for (const auto& element : elements(gridView))
ElementIterator end = gridView.template end<0>();
for (; it != end; ++it)
{ {
const Element& e = *it; const Geometry geometry = element.geometry();
const Geometry geometry = e.geometry(); const FiniteElement& fe = cache.get(element.type());
const FiniteElement& fe = cache.get(e.type());
int localSize = fe.localBasis().size(); int localSize = fe.size();
// get quadrature rule of appropiate order (P1/Q1) // get quadrature rule of appropiate order (P1/Q1)
int order = 0; int order = (element.type().isSimplex())
if (e.type().isSimplex()) ? 2*1
order = 2*1; : 2*dim;
else
order = 2*dim; const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
const QuadratureRule& quad = Dune::template QuadratureRules<double, dim>::rule(e.type(), order);
// store values of shape functions // store values of shape functions
std::vector<RangeType> values(localSize); std::vector<RangeType> values(localSize);
...@@ -287,7 +256,7 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f) ...@@ -287,7 +256,7 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
double z = quad[pt].weight() * integrationElement; double z = quad[pt].weight() * integrationElement;
for (int i = 0; i < localSize; ++i) for (int i = 0; i < localSize; ++i)
{ {
int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim); int iGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(i).subEntity(), dim);
r[iGlobal].axpy(values[i]*z, fAtPos); r[iGlobal].axpy(values[i]*z, fAtPos);
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment