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

Merge branch 'modernize-operatorassembler' into 'master'

Allow to set the partition type of the OperatorAssembler class

See merge request agnumpde/dune-fufem!48
parents 7e1e1b3f 3d001cb3
No related branches found
No related tags found
No related merge requests found
#ifndef OPERATOR_ASSEMBLER_HH #ifndef OPERATOR_ASSEMBLER_HH
#define OPERATOR_ASSEMBLER_HH #define OPERATOR_ASSEMBLER_HH
#include <dune/grid/common/rangegenerators.hh>
#include <dune/grid/common/partitionset.hh>
#include <dune/istl/matrix.hh> #include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
...@@ -8,8 +11,11 @@ ...@@ -8,8 +11,11 @@
#include "dune/fufem/functionspacebases/functionspacebasis.hh" #include "dune/fufem/functionspacebases/functionspacebasis.hh"
//! Generic global assembler for operators on a gridview /** \brief Generic global assembler for operators on a gridview
template <class TrialBasis, class AnsatzBasis> *
* \tparam PartitionSet For distributed grids: The subset of the grid elements to assemble on
*/
template <class TrialBasis, class AnsatzBasis, class PartitionSet=Dune::Partitions::All>
class OperatorAssembler class OperatorAssembler
{ {
private: private:
...@@ -67,24 +73,21 @@ class OperatorAssembler ...@@ -67,24 +73,21 @@ class OperatorAssembler
template <class LocalAssemblerType, bool lumping> template <class LocalAssemblerType, bool lumping>
void addIndicesStaticLumping(LocalAssemblerType& localAssembler, Dune::MatrixIndexSet& indices) const void addIndicesStaticLumping(LocalAssemblerType& localAssembler, Dune::MatrixIndexSet& indices) const
{ {
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename LocalAssemblerType::BoolMatrix BoolMatrix; typedef typename LocalAssemblerType::BoolMatrix BoolMatrix;
typedef typename TrialBasis::LinearCombination LinearCombination; typedef typename TrialBasis::LinearCombination LinearCombination;
ElementIterator it = tBasis_.getGridView().template begin<0>(); for (const auto& element : Dune::elements(tBasis_.getGridView(), PartitionSet()))
ElementIterator end = tBasis_.getGridView().template end<0>();
for (; it != end; ++it)
{ {
// get shape functions // get shape functions
const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(*it); const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(element);
const typename AnsatzBasis::LocalFiniteElement& aFE = aBasis_.getLocalFiniteElement(*it); const typename AnsatzBasis::LocalFiniteElement& aFE = aBasis_.getLocalFiniteElement(element);
BoolMatrix localIndices(tFE.localBasis().size(), aFE.localBasis().size()); BoolMatrix localIndices(tFE.localBasis().size(), aFE.localBasis().size());
localAssembler.indices(*it, localIndices, tFE, aFE); localAssembler.indices(element, localIndices, tFE, aFE);
for (size_t i=0; i<tFE.localBasis().size(); ++i) for (size_t i=0; i<tFE.localBasis().size(); ++i)
{ {
int rowIndex = tBasis_.index(*it, i); int rowIndex = tBasis_.index(element, i);
const LinearCombination& rowConstraints = tBasis_.constraints(rowIndex); const LinearCombination& rowConstraints = tBasis_.constraints(rowIndex);
bool rowIsConstrained = tBasis_.isConstrained(rowIndex); bool rowIsConstrained = tBasis_.isConstrained(rowIndex);
for (size_t j=0; j<aFE.localBasis().size(); ++j) for (size_t j=0; j<aFE.localBasis().size(); ++j)
...@@ -103,7 +106,7 @@ class OperatorAssembler ...@@ -103,7 +106,7 @@ class OperatorAssembler
} }
else else
{ {
int colIndex = aBasis_.index(*it, j); int colIndex = aBasis_.index(element, j);
const LinearCombination& colConstraints = aBasis_.constraints(colIndex); const LinearCombination& colConstraints = aBasis_.constraints(colIndex);
bool colIsConstrained = aBasis_.isConstrained(colIndex); bool colIsConstrained = aBasis_.isConstrained(colIndex);
if (rowIsConstrained and colIsConstrained) if (rowIsConstrained and colIsConstrained)
...@@ -137,24 +140,21 @@ class OperatorAssembler ...@@ -137,24 +140,21 @@ class OperatorAssembler
template <class LocalAssemblerType, class GlobalMatrixType, bool lumping> template <class LocalAssemblerType, class GlobalMatrixType, bool lumping>
void addEntriesStaticLumping(LocalAssemblerType& localAssembler, GlobalMatrixType& A) const void addEntriesStaticLumping(LocalAssemblerType& localAssembler, GlobalMatrixType& A) const
{ {
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename LocalAssemblerType::LocalMatrix LocalMatrix; typedef typename LocalAssemblerType::LocalMatrix LocalMatrix;
typedef typename TrialBasis::LinearCombination LinearCombination; typedef typename TrialBasis::LinearCombination LinearCombination;
ElementIterator it = tBasis_.getGridView().template begin<0>(); for (const auto& element : Dune::elements(tBasis_.getGridView(), PartitionSet()))
ElementIterator end = tBasis_.getGridView().template end<0>();
for (; it != end; ++it)
{ {
// get shape functions // get shape functions
const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(*it); const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(element);
const typename AnsatzBasis::LocalFiniteElement& aFE = aBasis_.getLocalFiniteElement(*it); const typename AnsatzBasis::LocalFiniteElement& aFE = aBasis_.getLocalFiniteElement(element);
LocalMatrix localA(tFE.localBasis().size(), aFE.localBasis().size()); LocalMatrix localA(tFE.localBasis().size(), aFE.localBasis().size());
localAssembler.assemble(*it, localA, tFE, aFE); localAssembler.assemble(element, localA, tFE, aFE);
for (size_t i=0; i<tFE.localBasis().size(); ++i) for (size_t i=0; i<tFE.localBasis().size(); ++i)
{ {
int rowIndex = tBasis_.index(*it, i); int rowIndex = tBasis_.index(element, i);
const LinearCombination& rowConstraints = tBasis_.constraints(rowIndex); const LinearCombination& rowConstraints = tBasis_.constraints(rowIndex);
bool rowIsConstrained = tBasis_.isConstrained(rowIndex); bool rowIsConstrained = tBasis_.isConstrained(rowIndex);
for (size_t j=0; j<aFE.localBasis().size(); ++j) for (size_t j=0; j<aFE.localBasis().size(); ++j)
...@@ -174,7 +174,7 @@ class OperatorAssembler ...@@ -174,7 +174,7 @@ class OperatorAssembler
} }
else else
{ {
int colIndex = aBasis_.index(*it, j); int colIndex = aBasis_.index(element, j);
const LinearCombination& colConstraints = aBasis_.constraints(colIndex); const LinearCombination& colConstraints = aBasis_.constraints(colIndex);
bool colIsConstrained = aBasis_.isConstrained(colIndex); bool colIsConstrained = aBasis_.isConstrained(colIndex);
if (rowIsConstrained and colIsConstrained) if (rowIsConstrained and colIsConstrained)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment