Skip to content
Snippets Groups Projects
Commit a892646e authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Some C++ modernization

parent c2eb7c82
No related branches found
No related tags found
No related merge requests found
...@@ -15,7 +15,6 @@ template <class Basis, class VectorType> ...@@ -15,7 +15,6 @@ template <class Basis, class VectorType>
class FEAssembler { class FEAssembler {
typedef typename Basis::GridView GridView; typedef typename Basis::GridView GridView;
typedef typename GridView::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator ElementIterator;
//! Dimension of the grid. //! Dimension of the grid.
enum { gridDim = GridView::dimension }; enum { gridDim = GridView::dimension };
...@@ -72,19 +71,16 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const ...@@ -72,19 +71,16 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
nb.resize(n, n); nb.resize(n, n);
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>(); for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition> (); {
const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(element);
for (; it!=endit; ++it) {
const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(*it);
for (size_t i=0; i<lfe.localBasis().size(); i++) { for (size_t i=0; i<lfe.localBasis().size(); i++) {
for (size_t j=0; j<lfe.localBasis().size(); j++) { for (size_t j=0; j<lfe.localBasis().size(); j++) {
int iIdx = basis_.index(*it,i); int iIdx = basis_.index(element,i);
int jIdx = basis_.index(*it,j); int jIdx = basis_.index(element,j);
nb.add(iIdx, jIdx); nb.add(iIdx, jIdx);
...@@ -116,33 +112,30 @@ assembleGradientAndHessian(const VectorType& sol, ...@@ -116,33 +112,30 @@ assembleGradientAndHessian(const VectorType& sol,
gradient.resize(sol.size()); gradient.resize(sol.size());
gradient = 0; gradient = 0;
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>(); for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition> (); {
const int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
for( ; it != endit; ++it ) {
const int numOfBaseFct = basis_.getLocalFiniteElement(*it).localBasis().size();
// Extract local solution // Extract local solution
VectorType localSolution(numOfBaseFct); VectorType localSolution(numOfBaseFct);
VectorType localPointLoads(numOfBaseFct); VectorType localPointLoads(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++) for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[basis_.index(*it,i)]; localSolution[i] = sol[basis_.index(element,i)];
VectorType localGradient(numOfBaseFct); VectorType localGradient(numOfBaseFct);
// setup local matrix and gradient // setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient); localStiffness_->assembleGradientAndHessian(element, basis_.getLocalFiniteElement(element), localSolution, localGradient);
// Add element matrix to global stiffness matrix // Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) { for(int i=0; i<numOfBaseFct; i++) {
int row = basis_.index(*it,i); int row = basis_.index(element,i);
for (int j=0; j<numOfBaseFct; j++ ) { for (int j=0; j<numOfBaseFct; j++ ) {
int col = basis_.index(*it,j); int col = basis_.index(element,j);
hessian[row][col] += localStiffness_->A_[i][j]; hessian[row][col] += localStiffness_->A_[i][j];
} }
...@@ -150,7 +143,7 @@ assembleGradientAndHessian(const VectorType& sol, ...@@ -150,7 +143,7 @@ assembleGradientAndHessian(const VectorType& sol,
// Add local gradient to global gradient // Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++) for (int i=0; i<numOfBaseFct; i++)
gradient[basis_.index(*it,i)] += localGradient[i]; gradient[basis_.index(element,i)] += localGradient[i];
} }
...@@ -166,21 +159,18 @@ computeEnergy(const VectorType& sol) const ...@@ -166,21 +159,18 @@ computeEnergy(const VectorType& sol) const
if (sol.size()!=basis_.size()) if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!"); DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.getGridView().template end<0,Dune::Interior_Partition>();
// Loop over all elements // Loop over all elements
for (; it!=endIt; ++it) { for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
{
// Number of degrees of freedom on this element // Number of degrees of freedom on this element
size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size(); size_t nDofs = basis_.getLocalFiniteElement(element).localBasis().size();
VectorType localSolution(nDofs); VectorType localSolution(nDofs);
for (size_t i=0; i<nDofs; i++) for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[basis_.index(*it,i)]; localSolution[i] = sol[basis_.index(element,i)];
energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution); energy += localStiffness_->energy(element, basis_.getLocalFiniteElement(element), localSolution);
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment