ellipticfemscheme.hh

#ifndef __ELLIPTIC_FEMSCHEME_HH__
#define __ELLIPTIC_FEMSCHEME_HH__
 
// iostream includes
#include <iostream>
 
// include grid part
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
 
// include discrete function space
#include <dune/fem/space/lagrange.hh>
 
// adaptation ...
#include <dune/fem/function/adaptivefunction.hh>
#include <dune/fem/space/common/adaptmanager.hh>
 
// include discrete function
#include <dune/fem/function/blockvectorfunction.hh>
 
// Newton's iteration
#include <dune/fem/solver/newtoninverseoperator.hh>
 
// lagrange interpolation
#include <dune/fem/space/common/interpolate.hh>
 
// include norms
#include <dune/fem/misc/l2norm.hh>
#include <dune/fem/misc/h1norm.hh>
 
// include parameter handling
#include <dune/fem/io/parameter.hh>
 
// include output
#include <dune/fem/io/file/dataoutput.hh>
 
// local includes
#include "../algorithms/femschemeinterface.hh"
#include "../models/modelinterface.hh"
#include "../operators/ellipticoperator.hh"
#include "../operators/functionals/l2innerproductfunctional.hh"
#include "../operators/functionals/l2boundaryfunctional.hh"
#include "../operators/functionals/functionalexpression.hh"
#include "../operators/constraints/dirichletconstraints.hh"
#include "../functions/gridfunctionexpression.hh"
#include "../functions/boundaryfunctionexpression.hh"
 
#include "../estimators/ellipticestimator.hh"
#include "../algorithms/marking.hh"
#include "../common/dataoutput.hh"
 
// Select an appropriate solver, depending on ModelType and solver-family (ISTL, PESC ...)
#include "../common/solverselector.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunction, class Model, class InitialGuess>
    class EllipticFemSchemeBase
      : public AdaptiveFemScheme
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
 
      typedef typename DiscreteFunctionType::GridType GridType;
 
      typedef typename DiscreteFunctionType::GridPartType GridPartType;
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      typedef Model ModelImplementationType;
      typedef ModelInterface<ModelImplementationType> ModelType;
 
      typedef typename ModelType::OperatorPartsType OperatorPartsType;
 
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
 
      typedef SolverSelector<DiscreteFunctionType, OperatorPartsType> SolverSelectorType;
      typedef typename SolverSelectorType::LinearOperatorType LinearOperatorType;
      typedef typename SolverSelectorType::LinearInverseOperatorType LinearInverseOperatorType;
 
      typedef
      Fem::NewtonInverseOperator<LinearOperatorType, LinearInverseOperatorType>
      NonLinearInverseOperatorType;
 
      typedef Fem::RestrictProlongDefault<DiscreteFunctionType> RestrictionProlongationType;
 
      typedef Fem::AdaptationManager<GridType, RestrictionProlongationType> AdaptationManagerType;
 
      typedef typename ModelType::DirichletBoundaryFunctionType DirichletBoundaryFunctionType;
 
      typedef typename ModelType::DirichletWeightFunctionType DirichletWeightFunctionType;
 
      typedef
      decltype(std::declval<DirichletBoundaryFunctionType>()
               /
               std::declval<typename DirichletWeightFunctionType::GridFunctionType>())
        EffectiveDirichletFunctionType;
 
      typedef typename ModelType::DirichletIndicatorType DirichletIndicatorType;
 
      typedef DirichletConstraints<LinearOperatorType, EffectiveDirichletFunctionType, DirichletIndicatorType> ConstraintsOperatorType;
 
      typedef DifferentiableEllipticOperator<LinearOperatorType, OperatorPartsType, ConstraintsOperatorType> OperatorType;
 
      typedef typename ModelType::BulkForcesFunctionType BulkForcesFunctionType;
      typedef
      L2InnerProductFunctional<DiscreteFunctionSpaceType, BulkForcesFunctionType>
      BulkForcesFunctional;
      typedef typename ModelType::NeumannBoundaryFunctionType BoundaryFluxFunctionType;
      typedef
      L2BoundaryFunctional<DiscreteFunctionSpaceType, BoundaryFluxFunctionType, typename ModelType::NeumannIndicatorType>
      BoundaryFluxFunctional;
 
      typedef EllipticEstimator<DiscreteFunctionSpaceType, ModelImplementationType> EstimatorType ;
 
      typedef MarkingStrategy<EstimatorType> MarkingStrategyType;
 
      typedef InitialGuess InitialGuessType;
 
      typedef
      typename GridFunctionConverter<InitialGuessType, GridPartType>::WrappedGridFunctionType
      InitialGuessFunctionType;
 
      typedef
      std::tuple<DiscreteFunctionType *, const InitialGuessFunctionType *>
      IOTupleType;
 
      typedef DataOutput<GridType, IOTupleType> DataOutputType;
 
     protected:
      EllipticFemSchemeBase(DiscreteFunctionType& solution,
                            const ModelType& model,
                            const InitialGuessType& initialGuess,
                            const std::string& name = "ellipt")
        : grid_(solution.space().gridPart().grid()),
          model_(model),
          name_(name),
          gridPart_(solution.space().gridPart()),
          discreteSpace_(solution.space()),
          solution_(solution),
          // estimator
          estimator_(model_(), discreteSpace_),
          markingStrategy_(estimator_),
          // restriction/prolongation operator
          restrictProlong_(0),
          // adaptation manager
          adaptationManager_(0),
          // create Dirichlet contraints
          boundaryValues_(model_().dirichletBoundaryFunction(gridPart_)),
          boundaryWeight_(model_().dirichletWeightFunction(gridPart_)),
          constraints_(discreteSpace_,
                       boundaryValues_ / boundaryWeight_.function(),
                       model_().dirichletIndicator()),
          // the elliptic operator (implicit)
          operator_(model_().operatorParts(), constraints_),
          // create linear operator (domainSpace,rangeSpace)
          linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
          solverEps_(Fem::Parameter::getValue<double>(name + ".solvereps", 1e-8)),
          sequence_(-1),
          // the right-hand-side
          bulkForces_(model_().bulkForcesFunction(gridPart_)),
          boundaryFlux_(model_().neumannBoundaryFunction(gridPart_)),
          // exact solution, if any
          initialGuess_(wrapToGridFunction("Initial Guess/Exact Solution", initialGuess, gridPart_, discreteSpace_.order()+1)),
          // io tuple
          ioTuple_(&solution_, &initialGuess_()),
          // DataOutputType
          dataOutput_(0)
      {}
 
#if 0
      // too complicated
      EllipticFemSchemeBase(const EllipticFemSchemeBase& other)
        : grid_(other.grid_),
          model_(other.model_),
          gridPart_(other.gridPart_),
          discreteSpace_(other.discreteSpace_),
          solution_(other.solution_),
          // estimator
          estimator_(other.estimator_),
          markingStrategy_(estimator_),
          // restriction/prolongation operator
          restrictProlong_(0),
          // adaptation manager
          adaptationManager_(0),
          // create Dirichlet contraints
          boundaryValues_(model_().dirichletBoundaryFunction(gridPart_)),
          boundaryWeight_(model_().dirichletWeightFunction(gridPart_)),
          constraints_(other.constraints_),
          // the elliptic operator (implicit)
          operator_(model_().operatorParts(), constraints_),
          // create linear operator (domainSpace,rangeSpace)
          linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
          solverEps_(other.solverEps_),
          sequence_(-1),
          // the right-hand-side
          bulkForces_(model_().bulkForcesFunction(gridPart_)),
          boundaryFlux_(model_().neumannBoundaryFunction(gridPart_)),
          // exact solution, if any
          initialGuess_(other.initialGuess_),
          // io tuple
          ioTuple_(&solution_, &initialGuess_),
          // DataOutputType
          dataOutput_(0)
      {}
#endif
 
     public:
 
      virtual ~EllipticFemSchemeBase() {
        if (dataOutput_) {
          delete dataOutput_;
        }
        if (adaptationManager_) {
          delete adaptationManager_;
        }
        if (restrictProlong_) {
          delete restrictProlong_;
        }
      }
 
      virtual void initialize()
      {
        if (std::is_same<
            InitialGuessFunctionType,
            ZeroGridFunctionTraits<FunctionSpaceType, GridPartType> >::value) {
          solution_.clear();
        } else {
          // apply natural interpolation
          interpolate(initialGuess_(), solution_);
        }
      }
 
      virtual void solve(bool forceMatrixAssembling = true)
      {
        // set Dirichlet constraints in solution
        constraints_.constrain(solution_);
 
        // right hand side, if present
        DiscreteFunctionType rhs("rhs", discreteSpace_);
 
        // If we have a RHS _function_ then also take that into
        // account. Note that it would be comparatively cheap to
        // omit the enclosing if-claus, because the default
        // configuration results into an efficient no-op. There is
        // also no point to optimize further into the case
        // Forces-but-no-Neumann or Neumann-but-no-forces for the
        // same reason.
        auto rhsFunctional =
          BulkForcesFunctional(discreteSpace_, bulkForces_)
          +
          BoundaryFluxFunctional(discreteSpace_, boundaryFlux_, model_().neumannIndicator())
          +
          model_().forcesFunctional(discreteSpace_);
 
        if (!isZero(rhsFunctional)) { // avoid copying thousands of zero values
          rhsFunctional.coefficients(rhs);
          constraints_.zeroConstrain(rhs);
        } else {
          rhs.clear();
        }
 
        if (OperatorPartsType::isLinear) {
          linearSolve(rhs, forceMatrixAssembling);
        } else {
          nonLinearSolve(rhs);
        }
      }
 
     protected:
 
      virtual void nonLinearSolve(DiscreteFunctionType& rhs)
      {
        assert(!OperatorPartsType::isLinear);
 
        NonLinearInverseOperatorType newton(operator_);
 
        newton(rhs, solution_);
 
        assert(newton.converged());
      }
 
      virtual void linearSolve(DiscreteFunctionType& rhs, bool forceMatrixAssembling)
      {
        assert(OperatorPartsType::isLinear);
 
        // if sequence number is outdated update linear operator
        if (sequence_ != discreteSpace_.sequence() || forceMatrixAssembling) {
          // assemble linear operator (i.e. setup matrix) The jacobian
          // incorporates the constraints defined by the constraint class.
          operator_.jacobian(solution_ , linearOperator_);
 
          // update sequence number
          sequence_ = discreteSpace_.sequence();
        }
 
        // inverse operator using linear operator
        LinearInverseOperatorType solver(linearOperator_, solverEps_, solverEps_);
 
        // Apply the affine linear operator to the start value. This
        // computes the initial residual. In the linear case, Lagrange space,
        // Dirichlet data g, this computes
        //
        // update_ = A u - g
        //
        // where it is allowed that A is affine-linear, without having
        // to apply a non-linear solver for this trivial non-linearity
        DiscreteFunctionType update("update", discreteSpace_);
        operator_(solution_, update);
        rhs -= update; // rhs = rhs - A u ...
 
        // Zero initial values for the update_ vector are ok, because
        // the information about the previous solution is already is
        // contained in residual_.
        update.clear();
 
        solver(rhs, update);
        solution_ += update; // ... so u = u + invA(rhs - Au)
      }
 
   public:
 
      virtual bool mark (const double tolerance)
      {
        return markingStrategy_.mark(tolerance);
      }
 
      virtual double estimate()
      {
        return estimator_.estimate(solution_);
      }
 
      virtual void adapt()
      {
        // there can only one adaptation manager per grid, and the
        // RestrictionProlongationType which determines which
        // functions are restricted/prolongated is built into the
        // AdaptationManagerType. In order to allow for override by
        // derived classes, allocation that adaptationManager_
        // dynamically (and thus not at all, if ThisType::adapt() is
        // never called.
 
        if (!adaptationManager_) {
          // restriction/prolongation operator
          restrictProlong_ = new RestrictionProlongationType(solution_);
 
          // adaptation manager
          adaptationManager_ = new AdaptationManagerType(grid_, *restrictProlong_);
        }
 
        // apply adaptation and load balancing
        adaptationManager_->adapt();
        if (adaptationManager_->loadBalance()) {
          // TODO: re-initialize stuff as needed.
        }
      }
 
      virtual int output()
      {
        if (!dataOutput_) {
          // NOTE: this should allocated dynamically, otherwise a
          // derived class has problems to define completely different
          // IO-schemes Also DataOutputType likes to already generate
          // some files during construction, so make sure this never happens
          dataOutput_ = new DataOutputType(grid_, ioTuple_, DataOutputParameters());
        }
 
        if (!dataOutput_->willWrite()) {
          return -1;
        }
 
        // write data
        dataOutput_->write();
 
        return dataOutput_->writeStep() - 1;
      }
 
      virtual double residual() const
      {
        // right hand side, if present
        DiscreteFunctionType rhs("rhs", discreteSpace_);
        auto rhsFunctional =
          BulkForcesFunctional(discreteSpace_, bulkForces_)
          +
          BoundaryFluxFunctional(discreteSpace_, boundaryFlux_, model_().neumannIndicator())
          +
          model_().forcesFunctional(discreteSpace_);
 
        if (!isZero(rhsFunctional)) { // avoid copying thousands of zero values
          rhsFunctional.coefficients(rhs);
          constraints_.zeroConstrain(rhs);
        } else {
          rhs.clear();
        }
        DiscreteFunctionType update("update", discreteSpace_);
        operator_(solution_, update);
        rhs -= update; // rhs = rhs - A u ...
        return std::sqrt(rhs.scalarProductDofs(rhs));
      }
 
      virtual double error() const
      {
        // can also use L2-norm
        typedef Fem::H1Norm< GridPartType > NormType;
 
        // the DiscreteFunctionSpaceAdapter sets the order per default
        // to 111 == \infty. We set therefore the order just high
        // enough that we see the asymptotic error. If the polynomial
        // order it D, then the error should decay with (h^D). In
        // principle it should thus suffice to use a quadrature of
        // degree D.
        NormType norm(gridPart_, 2*discreteSpace_.order());
        return norm.distance(initialGuess_(), solution_);
      }
 
      virtual size_t size() const
      {
        // NOTE: slaveDofs.size() is always one TOO LARGE.
        size_t numberOfDofs = discreteSpace_.blockMapper().size() - discreteSpace_.slaveDofs().size() + 1;
 
        numberOfDofs = grid_.comm().sum(numberOfDofs);
 
        return numberOfDofs;
      }
 
     protected:
      GridType&     grid_; // hierarchical grid
      ExpressionStorage<ModelImplementationType> model_;   // the mathematical model
 
      const std::string name_;
 
      GridPartType& gridPart_; // grid part(view), here the leaf grid the discrete space is build with
 
      const DiscreteFunctionSpaceType& discreteSpace_; // discrete function space
      DiscreteFunctionType& solution_;   // the unknown
 
      EstimatorType       estimator_;  // residual error estimator
      MarkingStrategyType markingStrategy_;
 
      RestrictionProlongationType *restrictProlong_ ; // local restriction/prolongation object
      AdaptationManagerType  *adaptationManager_ ;    // adaptation manager handling adaptation
 
      DirichletBoundaryFunctionType boundaryValues_;
      DirichletWeightFunctionType boundaryWeight_;
      ConstraintsOperatorType constraints_; // dirichlet boundary constraints
 
      OperatorType operator_; // the affine-linear operator.
 
      LinearOperatorType linearOperator_;  // the linear operator (i.e. jacobian of the operator)
 
      const double solverEps_ ; // eps for linear solver
      mutable int sequence_;            // sequence number
 
      BulkForcesFunctionType bulkForces_;
      BoundaryFluxFunctionType boundaryFlux_;
 
      ExpressionStorage<InitialGuessFunctionType> initialGuess_;
      IOTupleType ioTuple_; // tuple with pointers
      DataOutputType *dataOutput_; // data output class
    };
 
    template<class DiscreteFunction,
             class Model,
             class InitialGuess = ZeroGridFunction<typename Model::FunctionSpaceType,
                                                   typename Model::GridPartType> >
    class EllipticFemScheme;
 
    template<class DiscreteFunction, class Model, class InitialGuess>
    class EllipticFemScheme
      : public EllipticFemSchemeBase<DiscreteFunction, Model, InitialGuess>
    {
      typedef EllipticFemSchemeBase<DiscreteFunction, Model, InitialGuess> BaseType;
     public:
      typedef typename BaseType::DiscreteFunctionType DiscreteFunctionType;
      typedef typename BaseType::ModelType ModelType;
      typedef typename BaseType::InitialGuessType InitialGuessType;
 
      EllipticFemScheme(DiscreteFunctionType& solution,
                        const ModelType& model,
                        const InitialGuessType& initialGuess,
                        const std::string& name = "ellipt")
        : BaseType(solution, model, initialGuess, name)
      {}
    };
 
    template<class DiscreteFunction, class Model>
    class EllipticFemScheme<DiscreteFunction,
                            Model,
                            ZeroGridFunction<typename Model::FunctionSpaceType,
                                             typename Model::GridPartType> >
      : public EllipticFemSchemeBase<DiscreteFunction, Model,
                                     ZeroGridFunction<typename Model::FunctionSpaceType,
                                                      typename Model::GridPartType> >
    {
     public:
      typedef
      ZeroGridFunction<typename Model::FunctionSpaceType, typename Model::GridPartType>
      InitialGuessType;
     private:
      typedef EllipticFemSchemeBase<DiscreteFunction, Model, InitialGuessType> BaseType;
     public:
      typedef typename BaseType::DiscreteFunctionType DiscreteFunctionType;
      typedef typename BaseType::ModelType ModelType;
 
      EllipticFemScheme(DiscreteFunctionType& solution,
                        const ModelType& model,
                        const InitialGuessType& initialGuess,
                        const std::string& name = "ellipt")
        : BaseType(solution, model, initialGuess, name)
      {}
 
      EllipticFemScheme(DiscreteFunctionType& solution,
                        const ModelType& model,
                        const std::string& name = "ellipt")
        : BaseType(solution, model, InitialGuessType(solution.space().gridPart()), name)
      {}
    };
 
#if 0
    // too complicated, the scheme is rather not copy constructible ATM
 
    template<class DiscreteFunction, class Model, class InitialGuess>
    EllipticFemScheme<DiscreteFunction, Model, InitialGuess>
    ellipticFemScheme(DiscreteFunction& solution,
                      const Model&model,
                      const InitialGuess& initialGuess)
    {
      typedef EllipticFemScheme<DiscreteFunction, Model, InitialGuess> ReturnType;
 
      return ReturnType(solution, model, initialGuess);
    }
 
    template<class DiscreteFunction, class Model>
    EllipticFemScheme<DiscreteFunction,
                      Model,
                      ZeroGridFunction<typename Model::FunctionSpaceType,
                                       typename Model::GridPartType> >
    ellipticFemScheme(DiscreteFunction& solution, const Model& model)
    {
      typedef
        EllipticFemScheme<DiscreteFunction,
                          Model,
                          ZeroGridFunction<typename Model::FunctionSpaceType,
                                           typename Model::GridPartType> >
        ReturnType;
 
      return ReturnType(solution, model);
    }
#endif
 
 
  } // ACFem::
 
} // Dune::
 
#endif // __ELLIPTIC_FEMSCHEME_HH__