massoperator.hh

#ifndef MASSOPERATOR_HH
#define MASSOPERATOR_HH
 
#include <dune/common/dynvector.hh>
#include <dune/fem/common/bindguard.hh>
#include <dune/fem/function/common/localcontribution.hh>
#include <dune/fem/function/localfunction/const.hh>
#include <dune/fem/operator/common/differentiableoperator.hh>
#include <dune/fem/operator/common/localcontribution.hh>
#include <dune/fem/operator/common/operator.hh>
#include <dune/fem/operator/common/stencil.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/grid/common/rangegenerators.hh>
 
template <class DiscreteFunction, class LinearOperator>
class MassOperator : public LinearOperator {
  typedef MassOperator<DiscreteFunction, LinearOperator> ThisType;
  typedef LinearOperator BaseType;
 
  typedef DiscreteFunction DiscreteFunctionType;
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType
      DiscreteFunctionSpaceType;
 
  typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
  typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
 
  typedef Dune::Fem::CachingQuadrature<
      typename DiscreteFunctionSpaceType::GridPartType, 0>
      QuadratureType;
 
 public:
  explicit MassOperator(const DiscreteFunctionSpaceType &dfSpace)
      : BaseType("mass operator", dfSpace, dfSpace), dfSpace_(dfSpace) {
    assemble(*this);
  }
 
  template <class Function>
  void assembleRHS(const Function &u, DiscreteFunctionType &w) const;
  void assemble(LinearOperator &matrix) const;
 
 private:
  const DiscreteFunctionSpaceType &dfSpace_;
};
 
template <class DiscreteFunction, class LinearOperator>
class AffineMassOperator
    : public Dune::Fem::DifferentiableOperator<LinearOperator> {
 public:
  typedef DiscreteFunction DiscreteFunctionType;
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType
      DiscreteFunctionSpaceType;
  typedef LinearOperator JacobianOperatorType;
  template <class Function>
  explicit AffineMassOperator(const DiscreteFunctionSpaceType &dfSpace,
                              const Function &u)
      : matrix_(dfSpace), rhs_("rhs", dfSpace) {
    matrix_.assembleRHS(u, rhs_);
  }
  const DiscreteFunctionType rhs() const { return rhs_; }
  virtual void operator()(const DiscreteFunctionType &u,
                          DiscreteFunctionType &w) const {
    matrix_(u, w);
    w -= rhs_;
  }
  virtual void jacobian(const DiscreteFunctionType &u,
                        JacobianOperatorType &jOp) const {
    matrix_.assemble(jOp);
  }
 
 private:
  MassOperator<DiscreteFunction, LinearOperator> matrix_;
  DiscreteFunction rhs_;
};
 
template <class DiscreteFunction, class LinearOperator>
template <class Function>
void MassOperator<DiscreteFunction, LinearOperator>::assembleRHS(
    const Function &u, DiscreteFunctionType &w) const {
  // clear result
  w.clear();
 
  Dune::Fem::ConstLocalFunction<Function> uLocal(u);
  Dune::Fem::AddLocalContribution<DiscreteFunctionType> wLocal(w);
 
  // run over entities
  for (const auto &entity :
       elements(dfSpace_.gridPart(), Dune::Partitions::interiorBorder)) {
    const auto geometry = entity.geometry();
 
    auto guard = bindGuard(std::tie(uLocal, wLocal), entity);
 
    // run over quadrature points
    for (const auto qp :
         QuadratureType(entity, uLocal.order() + wLocal.order())) {
      // evaluate u
      RangeType uValue = uLocal.evaluate(qp);
 
      // apply quadrature weight
      uValue *= qp.weight() * geometry.integrationElement(qp.position());
 
      // add to local function
      wLocal.axpy(qp, uValue);
    }
  }
}
 
template <class DiscreteFunction, class LinearOperator>
void MassOperator<DiscreteFunction, LinearOperator>::assemble(
    LinearOperator &matrix) const {
  matrix.reserve(Dune::Fem::DiagonalStencil<DiscreteFunctionSpaceType,
                                            DiscreteFunctionSpaceType>(
      dfSpace_, dfSpace_));
  matrix.clear();
 
  Dune::Fem::AddLocalContribution<LinearOperator> localMatrix(matrix);
  Dune::DynamicVector<typename DiscreteFunctionSpaceType::RangeType> values;
 
  // run over entities
  for (const auto &entity :
       elements(dfSpace_.gridPart(), Dune::Partitions::interiorBorder)) {
    const auto geometry = entity.geometry();
 
    auto guard = bindGuard(localMatrix, entity, entity);
 
    const auto &basis = localMatrix.domainBasisFunctionSet();
    const unsigned int numBasisFunctions = basis.size();
    values.resize(numBasisFunctions);
 
    // run over quadrature points
    for (const auto qp : QuadratureType(entity, localMatrix.order())) {
      // evaluate base functions
      basis.evaluateAll(qp, values);
 
      // apply quadrature weight
      const auto weight =
          qp.weight() * geometry.integrationElement(qp.position());
      std::for_each(values.begin(), values.end(),
                    [weight](auto &a) { a *= weight; });
 
      // update system matrix
      localMatrix.axpy(qp, values);
    }
  }
}
 
#endif  // #ifndef MASSOPERATOR_HH