feassembler.hh

#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 (size_t i=0; i<size; i++)
      localConfiguration[i] = configurationBackend[ localView.index(i) ];

    // setup local matrix and gradient
    localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian);

    for (size_t i=0; i<size; i++)
    {
      auto row = localView.index(i);
      gradientBackend[row] += localGradient[i];

      for (size_t 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:
    [[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