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

Merge branch 'port-to-newer-fufem' into 'master' 

Port TrustRegionSolver to current dune-fufem

See merge request ag-sander/dune/dune-elasticity!78
parents a684da1e fd5f3ff7
No related branches found
No related tags found
No related merge requests found
Pipeline #11202 failed
Stage: .pre
Stage: test
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
+ 3
− 0
View file @ f750db35
# 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
......
dune/elasticity/assemblers/feassembler.hh
+ 0
− 184
View file @ f750db35
......@@ -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
dune/elasticity/common/trustregionsolver.cc
+ 22
− 111
View file @ f750db35
......@@ -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_,
......
dune/elasticity/common/trustregionsolver.hh
+ 4
− 10
View file @ f750db35
......@@ -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_;
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment