diff --git a/CHANGELOG.md b/CHANGELOG.md
index ed9fe0f88cb8922604e52e23f45f53e6133abfb8..fc8aefdb120826cf92dd426b05c2ca052e7e2bc5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
# Master (will become release 2.10)
+- The class `Dune::FEAssembler` has been removed -- please use `Dune::Elasticity::FEAssembler`
+ from now on. The removed code still used the old `dune-fufem` function space basis interface.
+
# Release 2.9
diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh
index 8b922c6fc0f7eb43a6264899396d1b5734a18898..e7f062db00cf3573413aec7dc37092c38d2e8b18 100644
--- a/dune/elasticity/assemblers/feassembler.hh
+++ b/dune/elasticity/assemblers/feassembler.hh
@@ -218,188 +218,4 @@ computeEnergy(const VectorType& configuration) const
} // namespace Dune::Elasticity
-
-
-namespace Dune
-{
-
-/** \brief A global FE assembler for variational problems (old fufem bases version)
- */
-template <class Basis, class VectorType >
-class
-FEAssembler
-{
-
- typedef typename Basis::GridView GridView;
-
- //! Dimension of the grid.
- enum { gridDim = GridView::dimension };
-
- //! Dimension of a tangent space
- enum { blocksize = VectorType::value_type::dimension };
-
- //!
- typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
-public:
- [[deprecated("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.")]]
- const Basis basis_;
-
- /** \brief Partition type on which to assemble
- *
- * We want a fully nonoverlapping partition, and therefore need to skip ghost elements.
- */
- static constexpr auto partitionType = Dune::Partitions::interiorBorder;
-
-protected:
-
- LocalFEStiffness<GridView,
- typename Basis::LocalFiniteElement,
- VectorType>* localStiffness_;
-
-public:
-
- /** \brief Constructor for a given grid */
- FEAssembler(const Basis& basis,
- LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType>* localStiffness)
- : basis_(basis),
- localStiffness_(localStiffness)
- {}
-
- /** \brief Assemble the tangent stiffness matrix and the functional gradient together
- *
- * This may be more efficient than computing them separately
- */
- virtual void assembleGradientAndHessian(const VectorType& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
- Dune::BCRSMatrix<MatrixBlock>& hessian,
- bool computeOccupationPattern=true) const;
-
- /** \brief Compute the energy of a deformation state */
- virtual double computeEnergy(const VectorType& sol) const;
-
- //protected:
- void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
-
-}; // end class
-
-
-
-template <class Basis, class TargetSpace>
-void FEAssembler<Basis,TargetSpace>::
-getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
-{
- int n = basis_.size();
-
- nb.resize(n, n);
-
- for (const auto& element : elements(basis_.getGridView(), partitionType))
- {
- const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(element);
-
- for (size_t i=0; i<lfe.localBasis().size(); i++) {
-
- for (size_t j=0; j<lfe.localBasis().size(); j++) {
-
- int iIdx = basis_.index(element,i);
- int jIdx = basis_.index(element,j);
-
- nb.add(iIdx, jIdx);
-
- }
-
- }
-
- }
-
-}
-
-template <class Basis, class VectorType>
-void FEAssembler<Basis,VectorType>::
-assembleGradientAndHessian(const VectorType& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
- Dune::BCRSMatrix<MatrixBlock>& hessian,
- bool computeOccupationPattern) const
-{
- if (computeOccupationPattern) {
-
- Dune::MatrixIndexSet neighborsPerVertex;
- getNeighborsPerVertex(neighborsPerVertex);
- neighborsPerVertex.exportIdx(hessian);
-
- }
-
- hessian = 0;
-
- gradient.resize(sol.size());
- gradient = 0;
-
- for (const auto& element : elements(basis_.getGridView(), partitionType))
- {
- const int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
-
- // Extract local solution
- VectorType localSolution(numOfBaseFct);
- VectorType localPointLoads(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[basis_.index(element,i)];
-
- VectorType localGradient(numOfBaseFct);
-
- // setup local matrix and gradient
- localStiffness_->assembleGradientAndHessian(element, basis_.getLocalFiniteElement(element), localSolution, localGradient);
-
- // Add element matrix to global stiffness matrix
- for(int i=0; i<numOfBaseFct; i++) {
-
- int row = basis_.index(element,i);
-
- for (int j=0; j<numOfBaseFct; j++ ) {
-
- int col = basis_.index(element,j);
- hessian[row][col] += localStiffness_->A_[i][j];
-
- }
- }
-
- // Add local gradient to global gradient
- for (int i=0; i<numOfBaseFct; i++)
- gradient[basis_.index(element,i)] += localGradient[i];
-
- }
-
-}
-
-
-template <class Basis, class VectorType>
-double FEAssembler<Basis, VectorType>::
-computeEnergy(const VectorType& sol) const
-{
- double energy = 0;
-
- if (sol.size()!=basis_.size())
- DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
-
- // Loop over all elements
- for (const auto& element : elements(basis_.getGridView(), partitionType))
- {
- // Number of degrees of freedom on this element
- size_t nDofs = basis_.getLocalFiniteElement(element).localBasis().size();
-
- VectorType localSolution(nDofs);
-
- for (size_t i=0; i<nDofs; i++)
- localSolution[i] = sol[basis_.index(element,i)];
-
- energy += localStiffness_->energy(element, basis_.getLocalFiniteElement(element), localSolution);
-
- }
-
- return energy;
-
-}
-
-} // namespace Dune
-
#endif
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 7e9fb980840c64a4c1959603c0b92ea301357b47..ef683c5c07cc49f5d308c515fc89d8c902de15af 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -6,7 +6,6 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
@@ -20,8 +19,8 @@
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/twonorm.hh>
-#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/elasticity/common/maxnormtrustregion.hh>
@@ -29,10 +28,7 @@
template <class BasisType, class VectorType>
void TrustRegionSolver<BasisType,VectorType>::
setup(const typename BasisType::GridView::Grid& grid,
- std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
- Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
- const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
- const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler,
+ const Dune::Elasticity::FEAssembler<BasisType,VectorType>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -71,10 +67,7 @@ setup(const typename BasisType::GridView::Grid& grid,
BasisType feBasis(grid.levelGridView(grid.maxLevel()));
- std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
- Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
- BasisType,
- DuneFunctionsBasis<BasisType> > basis(feBasis);
+ BasisType basis(feBasis);
innerMultigridStep_ = std::make_unique<Dune::ParMG::Multigrid<VectorType> >();
innerMultigridStep_->mu_ = mu;
@@ -85,10 +78,7 @@ setup(const typename BasisType::GridView::Grid& grid,
std::cout << "Parallel multigrid setup took " << setupTimer.stop() << " seconds." << std::endl;
#else
- std::conditional_t< // do we have a dune-functions basis?
- Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
- BasisType,
- DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
+ BasisType basis(assembler_->basis_);
// ////////////////////////////////
// Create a multigrid solver
// ////////////////////////////////
@@ -131,123 +121,44 @@ setup(const typename BasisType::GridView::Grid& grid,
// //////////////////////////////////////////////////////////////////////////////////////
- using ScalarMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
- ScalarMatrixType localA;
+ MatrixType localA;
- std::conditional_t< // do we have a dune-functions basis?
- Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
- Dune::Fufem::DuneFunctionsOperatorAssembler<BasisType,BasisType>,
- OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> > operatorAssembler(basis,basis);
+ Dune::Fufem::DuneFunctionsOperatorAssembler operatorAssembler(basis,basis);
- if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
- {
- using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
- // construct a Fufem Basis Assembler
- auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
- // transform it to a dune-functions assembler
- auto localAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
- // the fufem basis assembler assumes a scalar basis but a blocked element matrix!
- // we can use ScaledIdentityMatrix here
- Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
- laplaceStiffness.assemble(element,
- localBlockedMatrix,
- trialLocalView.tree().child(0).finiteElement(),
- ansatzLocalView.tree().child(0).finiteElement());
- // convert back to flat matrix (we need only the diagonal entries in the blocks)
- localMatrix = 0;
- for(size_t i=0; i<localBlockedMatrix.N(); i++)
- for(size_t j=0; j<localBlockedMatrix.M(); j++)
- Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
- };
-
-
- auto matrixBackend = Dune::Fufem::istlMatrixBackend(localA);
- auto patternBuilder = matrixBackend.patternBuilder();
-
- operatorAssembler.assembleBulkPattern(patternBuilder);
- patternBuilder.setupMatrix();
-
- operatorAssembler.assembleBulkEntries(matrixBackend, localAssembler);
- }
- else
- {
- LaplaceAssembler<GridType,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> laplaceStiffness;
+ Dune::Fufem::LaplaceAssembler laplaceStiffness;
- operatorAssembler.assemble(laplaceStiffness, localA);
- }
+ operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localA), laplaceStiffness);
- ScalarMatrixType* A = new ScalarMatrixType(localA);
+ MatrixType* A = new MatrixType(localA);
- h1SemiNorm_ = std::make_shared< H1SemiNorm<CorrectionType> >(*A);
+ h1SemiNorm_ = std::make_shared<EnergyNorm<MatrixType, CorrectionType> >(*A);
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a mass matrix to create a norm that's equivalent to the L2-norm
// This will be used to monitor the gradient
// //////////////////////////////////////////////////////////////////////////////////////
- ScalarMatrixType localMassMatrix;
+ MatrixType localMassMatrix;
- if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
- {
- // construct a Fufem Basis Assembler
- using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
- auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
- // transform it to a dune-functions assembler
- auto localMassAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
- // the fufem basis assembler assumes a scalar basis but a blocked element matrix!
- // we can use ScaledIdentityMatrix here
- Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
- massStiffness.assemble(element,
- localBlockedMatrix,
- trialLocalView.tree().child(0).finiteElement(),
- ansatzLocalView.tree().child(0).finiteElement());
- // convert back to flat matrix (we need only the diagonal entries in the blocks)
- localMatrix = 0;
- for(size_t i=0; i<localBlockedMatrix.N(); i++)
- for(size_t j=0; j<localBlockedMatrix.M(); j++)
- Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
- };
-
- auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
- auto massPatternBuilder = massMatrixBackend.patternBuilder();
-
- operatorAssembler.assembleBulkPattern(massPatternBuilder);
- massPatternBuilder.setupMatrix();
-
- operatorAssembler.assembleBulkEntries(massMatrixBackend, localMassAssembler);
- }
- else
- {
- MassAssembler<GridType,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
- typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> massStiffness;
+ Dune::Fufem::MassAssembler massStiffness;
- operatorAssembler.assemble(massStiffness, localMassMatrix);
- }
+ operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localMassMatrix), massStiffness);
+
+ MatrixType* massMatrix = new MatrixType(localMassMatrix);
- ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
- l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
+ l2Norm_ = std::make_shared<EnergyNorm<MatrixType, CorrectionType> >(*massMatrix);
// ////////////////////////////////////////////////////////////
// Create Hessian matrix and its occupation structure
// ////////////////////////////////////////////////////////////
hessianMatrix_ = std::make_shared<MatrixType>();
- if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
- {
- auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
- auto hessianPatternBuilder = hessianBackend.patternBuilder();
- operatorAssembler.assembleBulkPattern(hessianPatternBuilder);
- hessianPatternBuilder.setupMatrix();
- }
- else
- {
- Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
- assembler_->getNeighborsPerVertex(indices);
- indices.exportIdx(*hessianMatrix_);
- }
+
+ auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
+ auto hessianPatternBuilder = hessianBackend.patternBuilder();
+ operatorAssembler.assembleBulkPattern(hessianPatternBuilder);
+ hessianPatternBuilder.setupMatrix();
+
#if !HAVE_DUNE_PARMG
innerSolver_ = std::make_shared<LoopSolver<CorrectionType> >(mmgStep,
innerIterations_,
diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh
index ca65bde1b302f6d6942bdcaadbc7e286494329ae..fcd49fde53ffc8f71aadc787e16249e111e691ad 100644
--- a/dune/elasticity/common/trustregionsolver.hh
+++ b/dune/elasticity/common/trustregionsolver.hh
@@ -77,10 +77,7 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
- std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
- Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
- const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
- const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>*> assembler,
+ const Dune::Elasticity::FEAssembler<BasisType,VectorType>* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -146,10 +143,7 @@ protected:
std::shared_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
- std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
- Dune::models< Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
- const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
- const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler_;
+ const Dune::Elasticity::FEAssembler<BasisType,VectorType>* assembler_;
#if HAVE_DUNE_PARMG
/** \brief The solver for the inner problem */
@@ -170,10 +164,10 @@ protected:
std::shared_ptr<const Dune::BitSetVector<blocksize> > ignoreNodes_;
/** \brief The norm used to measure multigrid convergence */
- std::shared_ptr< H1SemiNorm<CorrectionType> > h1SemiNorm_;
+ std::shared_ptr<EnergyNorm<MatrixType,CorrectionType> > h1SemiNorm_;
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */
- std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
+ std::shared_ptr<EnergyNorm<MatrixType,CorrectionType> > l2Norm_;
};