Select Git revision
lisa_julia.nebel_at_tu-dresden.de authored
Add missing include (grid/common/partitionset.hh) in feassembler.hh, as the FEAssembler uses Dune::Partitions::interiorBorder
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
feassembler.hh 10.34 KiB
#ifndef DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
#define DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
#include <dune/common/fmatrix.hh>
#include <dune/fufem/assemblers/istlbackend.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/grid/common/partitionset.hh>
#include "localfestiffness.hh"
namespace Dune::Elasticity {
/** \brief A global FE assembler for variational problems for dune-functions global bases
*/
template <class Basis, class VectorType>
class FEAssembler {
using LocalView = typename Basis::LocalView;
using field_type = typename VectorType::field_type;
//! Dimension of the grid.
enum { gridDim = LocalView::GridView::dimension };
public:
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:
std::shared_ptr<LocalFEStiffness<LocalView,field_type>> localStiffness_;
public:
/** \brief Constructor for a given basis and local assembler */
FEAssembler(const Basis& basis,
std::shared_ptr<LocalFEStiffness<LocalView, field_type>> 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
*/
template<class MatrixType>
void assembleGradientAndHessian(const VectorType& configuration,
VectorType& gradient,
MatrixType& hessian,
bool computeOccupationPattern=true) const;
/** \brief Compute the energy of a deformation state */
auto computeEnergy(const VectorType& configuration) const;
}; // end class
template <class Basis, class VectorType>
template <class MatrixType>
void FEAssembler<Basis,VectorType>::
assembleGradientAndHessian(const VectorType& configuration,
VectorType& gradient,
MatrixType& hessian,
bool computeOccupationPattern) const
{
// create backends for multi index access
auto hessianBackend = Dune::Fufem::ISTLMatrixBackend(hessian);
auto configurationBackend = Dune::Functions::istlVectorBackend(configuration);
auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
if (computeOccupationPattern)
{
auto patternBuilder = hessianBackend.patternBuilder();
patternBuilder.resize(basis_,basis_);
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView(), partitionType))
{
localView.bind(element);
for (size_t i=0; i<localView.size(); i++)
for (size_t j=0; j<localView.size(); j++)
patternBuilder.insertEntry( localView.index(i), localView.index(j) );
}
patternBuilder.setupMatrix();
}
gradientBackend.resize(basis_);
hessian = 0;
gradient = 0;
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView(), partitionType))
{
localView.bind(element);
const auto size = localView.size();
// Extract local configuration
std::vector<field_type> localConfiguration(size);
std::vector<field_type> localGradient(size);
Matrix<field_type> localHessian;
localHessian.setSize(size, size);
for (int i=0; i<size; i++)
localConfiguration[i] = configurationBackend[ localView.index(i) ];
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian);
for(int i=0; i<size; i++)
{
auto row = localView.index(i);
gradientBackend[row] += localGradient[i];
for (int j=0; j<size; j++ )
{
auto col = localView.index(j);
// fufem ISTLBackend uses operator() for entry access
hessianBackend(row,col) += localHessian[i][j];
}
}
}}
template <class Basis, class VectorType>
auto FEAssembler<Basis,VectorType>::
computeEnergy(const VectorType& configuration) const
{
double energy = 0;
if (configuration.dim()!=basis_.dimension())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
auto localView = basis_.localView();
// fufem multi index access
auto configurationBackend = Functions::istlVectorBackend(configuration);
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), partitionType))
{
localView.bind(element);
const auto size = localView.size();
std::vector<field_type> localConfiguration(size);
for (size_t i=0; i<size; i++)
localConfiguration[i] = configurationBackend[ localView.index(i) ];
energy += localStiffness_->energy(localView, localConfiguration);
}
return energy;
}
} // 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:
const Basis
DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.") 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