splittingfemscheme.hh

#ifndef __DUNE_ACFEM_SPLITTINGFEMSCHEME_HH__
#define __DUNE_ACFEM_SPLITTINGFEMSCHEME_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/operator/lagrangeinterpolation.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 vtk output and checkpointing
#include <dune/fem/io/file/dataoutput.hh>
#include <dune/fem/io/file/datawriter.hh>
 
// blah blah blah if only Dune were as hip as it ought to be
 
// local includes
#include "../algorithms/femschemeinterface.hh"
#include "../models/modelinterface.hh"
#include "../functions/gridfunctionwrapper.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 "../operators/constraints/emptyblockconstraints.hh"
#include "../estimators/parabolicestimator.hh"
#include "../estimators/trueerrorestimator.hh"
#include "../algorithms/marking.hh"
#include "../models/timeview.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 ImplicitModel,
             class ExplicitDataFunction,
             class ExplicitModel,
             class InitialGuessFunction>
    class SplittingFemSchemeBase
      : public BasicFemScheme
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
 
      typedef typename DiscreteFunctionType::GridType GridType;
 
      typedef typename DiscreteFunctionType::GridPartType GridPartType;
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      typedef ImplicitModel ImplicitModelType;
      typedef ModelConstituents<ImplicitModelType> ImplicitModelConstituents;
      typedef ExplicitModel ExplicitModelType;
      typedef ModelConstituents<ExplicitModelType> ExplicitModelConstituents;
 
      typedef
      decltype(std::declval<ImplicitModelType>() + std::declval<ExplicitModelType>())
        SumModelType;
      typedef ModelConstituents<SumModelType> SumModelConstituents;
 
      typedef ExplicitDataFunction ExplicitDataType;
      typedef
      typename GridFunctionConverter<ExplicitDataType, GridPartType>::WrappedGridFunctionType
      ExplicitDataFunctionType;
 
      typedef InitialGuessFunction InitialGuessType;
      typedef
      typename GridFunctionConverter<InitialGuessType, GridPartType>::WrappedGridFunctionType
      InitialGuessFunctionType;
 
      typedef typename ImplicitModelType::FunctionSpaceType FunctionSpaceType;
 
      // or: (must coincide)
      //typedef typenema DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
 
      // choose type of discrete function, Matrix implementation and solver implementation
      typedef SolverSelector<DiscreteFunctionType, ImplicitModelType> SolverSelectorType;
      typedef typename SolverSelectorType::LinearOperatorType LinearOperatorType;
      typedef typename SolverSelectorType::LinearInverseOperatorType LinearInverseOperatorType;
 
      typedef
      Fem::NewtonInverseOperator<LinearOperatorType, LinearInverseOperatorType>
      NonLinearInverseOperatorType;
 
      typedef Dune::Fem::RestrictProlongDefault<DiscreteFunctionType> RestrictionProlongationType;
 
      typedef Dune::Fem::AdaptationManager<GridType, RestrictionProlongationType> AdaptationManagerType;
 
      typedef typename ImplicitModelType::DirichletIndicatorType DirichletIndicatorType;
 
      typedef typename ImplicitModelType::DirichletBoundaryFunctionType DirichletBoundaryFunctionType;
 
      typedef typename ImplicitModelType::DirichletWeightFunctionType DirichletWeightFunctionType;
 
      typedef
      decltype(std::declval<DirichletBoundaryFunctionType>()
               /
               std::declval<typename DirichletWeightFunctionType::GridFunctionType>())
        EffectiveDirichletFunctionType;
 
      typedef DirichletConstraints<LinearOperatorType, EffectiveDirichletFunctionType, DirichletIndicatorType> ConstraintsOperatorType;
 
      typedef DifferentiableEmptyBlockConstraints<LinearOperatorType>  EmptyConstraintsType;
 
      typedef DifferentiableEllipticOperator<LinearOperatorType, ImplicitModelType, ConstraintsOperatorType> ImplicitOperatorType;
 
      typedef EllipticOperator<ExplicitModelType,
                               ExplicitDataFunctionType, DiscreteFunctionType,
                               EmptyConstraintsType> ExplicitOperatorType;
 
      // Isn't ther a more simpler way for doing this?
      typedef typename SumModelType::BulkForcesFunctionType BulkForcesFunctionType;
      typedef
      L2InnerProductFunctional<DiscreteFunctionSpaceType, BulkForcesFunctionType>
      BulkForcesFunctional;
      typedef typename SumModelType::NeumannBoundaryFunctionType BoundaryFluxFunctionType;
      typedef
      L2BoundaryFunctional<DiscreteFunctionSpaceType,
                           BoundaryFluxFunctionType,
                           typename ImplicitModelType::NeumannIndicatorType>
      BoundaryFluxFunctional;
 
      typedef InitialGuessFunctionType ExactSolutionFunctionType; // slight abuse
 
      typedef std::tuple<DiscreteFunctionType *, const ExactSolutionFunctionType *> IOTupleType;
 
      typedef DataOutput<GridType, IOTupleType> DataOutputType;
 
      typedef Dune::Fem::CheckPointer<GridType>   CheckPointerType;
 
      SplittingFemSchemeBase(DiscreteFunctionType& solution,
                             const ModelInterface<ImplicitModelType>& implicitModel,
                             const ExplicitDataType& explicitData,
                             const ModelInterface<ExplicitModelType>& explicitModel,
                             const InitialGuessType& initialGuess,
                             const std::string name = "femscheme")
        : solution_(solution)
        , implicitModel_(asImp(implicitModel))
        , explicitData_(explicitData)
        , explicitModel_(asImp(explicitModel))
        , initialGuess_(initialGuess)
        , name_(name)
        , grid_(solution.space().gridPart().grid())
        , gridPart_(solution.space().gridPart())
        , discreteSpace_(solution.space())
        , residual_("residual", discreteSpace_)
        , update_("update", discreteSpace_)
        // create Dirichlet contraints
        , boundaryValues_(implicitModel_().dirichletBoundaryFunction(gridPart_))
        , boundaryWeight_(implicitModel_().dirichletWeightFunction(gridPart_))
        , dirichletConstraints_(discreteSpace_,
                                boundaryValues_ / boundaryWeight_.function(),
                                implicitModel_().dirichletIndicator())
        // the elliptic operator (implicit)
        , implicitOperator_(implicitModel_().operatorParts(), dirichletConstraints_)
        // the elliptic operator (explicit)
        , explicitOperator_(explicitModel_().operatorParts())
        // create linear operator (domainSpace,rangeSpace)
        , linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_)
        , solverEps_(Dune::Fem::Parameter::getValue<double>(name_ + ".solvereps", 1e-8))
        // the right-hand-side
        , sumModel_(implicitModel_() + explicitModel_())
        , bulkForces_(sumModel_.bulkForcesFunction(gridPart_))
        , boundaryFlux_(sumModel_.neumannBoundaryFunction(gridPart_))
        // other stuff
        , sequence_(-1)
        , ioTuple_(&solution_, &initialGuess_())
        , dataOutput_(0)
      {}
 
      virtual ~SplittingFemSchemeBase() {
        if (dataOutput_) {
          delete dataOutput_;
        }
      }
 
      virtual void initialize()
      {
        // select Lagrange interpolation
        Fem::LagrangeInterpolation<InitialGuessFunctionType, DiscreteFunctionType> interpolation;
 
        if (std::is_same<
            InitialGuessFunctionType,
            ZeroGridFunctionTraits<FunctionSpaceType, GridPartType> >::value) {
          solution_.clear();
        } else {
          // apply lagrange interpolation
          interpolation(initialGuess_(), solution_);
        }
      }
 
      virtual void solve(bool forceMatrixAssembling = true)
      {
        // 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_, implicitModel_().neumannIndicator())
          +
          sumModel_.forcesFunctional(discreteSpace_);
 
        if (!isZero(rhsFunctional)) { // avoid copying thousands of zero values
          rhsFunctional.coefficients(residual_);
        } else {
          residual_.clear();
        }
 
        // Compute the contribution from the time derivative
        explicitOperator_(explicitData_(), update_);
        residual_ -= update_;
 
        if (ImplicitModelConstituents::hasDirichletBoundary) {
          // set Dirichlet constraints, new time
          dirichletConstraints_.rebuildValues(); // force update
          dirichletConstraints_.constrain(solution_);
 
          // do no spoil the new Dirichlet-values
          dirichletConstraints_.zeroConstrain(residual_);
        }
 
        if (ImplicitModelType::isLinear) {
          linearSolve(residual_, forceMatrixAssembling);    // sove Lu = res
        } else {
          nonLinearSolve(residual_); // solve Lu = res
        }
      }
 
     protected:
 
      virtual void nonLinearSolve(DiscreteFunctionType& rhs)
      {
        assert(!ImplicitModelType::isLinear);
 
        NonLinearInverseOperatorType newton(implicitOperator_);
 
        newton(rhs, solution_);
 
        assert(newton.converged());
      }
 
      virtual void linearSolve(DiscreteFunctionType& rhs, bool forceMatrixAssembling)
      {
        assert(ImplicitModelType::isLinear);
 
        if (sequence_ != discreteSpace_.sequence() || forceMatrixAssembling) {
          // assemble linear operator (i.e. setup matrix) The jacobian
          // incorporates the constraints defined by the constraint class.
          implicitOperator_.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
        implicitOperator_(solution_, update_);
        rhs -= update_; // rhs = rhs - A u ...
 
        // solve system. Zero initial values for the update_ vector are
        //  ok, because the information about the previous solution already
        // is contained in residual_.
        update_.clear();
        solver(rhs, update_);
 
        // subtract the update, difference should be the solution.
        solution_ += update_; // ... so u = u + invA(rhs - Au)
      }
 
     public:
 
      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;
        }
 
        dataOutput_->write();
 
        return dataOutput_->writeStep() - 1;
      }
 
      virtual double residual() const
      {
        assert(0);
        return -1;
      }
 
      virtual double error() const
      {
        // can also use H1-norm
        typedef Dune::Fem::L2Norm<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 is 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()+2);
        return norm.distance(initialGuess_(), solution_);
      }
 
      virtual size_t size() const
      {
        typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
        typedef Dune::Fem::SlaveDofs<DiscreteFunctionSpaceType, BlockMapperType> SlaveDofsType;
        typedef typename SlaveDofsType::SingletonKey SlaveDofsKeyType;
        typedef Dune::Fem::SingletonList<SlaveDofsKeyType, SlaveDofsType> SlaveDofsProviderType;
 
        SlaveDofsKeyType key(discreteSpace_, discreteSpace_.blockMapper());
        SlaveDofsType& slaveDofs(SlaveDofsProviderType::getObject(key));
 
        slaveDofs.rebuild();
 
        // slaveDofs.size() is off by one
        size_t numberOfDofs = discreteSpace_.blockMapper().size() - slaveDofs.size() + 1;
 
        numberOfDofs = grid_.comm().sum(numberOfDofs);
 
        numberOfDofs *= DiscreteFunctionSpaceType::localBlockSize;
 
        return numberOfDofs;
      }
 
     protected:
      DiscreteFunctionType& solution_;   // the unknown
      ExpressionStorage<ImplicitModelType> implicitModel_;
      ExpressionStorage<ExplicitDataFunctionType> explicitData_;  // this is already the GRID function
      ExpressionStorage<ExplicitModelType> explicitModel_; // old time-step
      ExpressionStorage<InitialGuessFunctionType> initialGuess_;  // this is already the GRID function
      const std::string name_;
 
      GridType&     grid_; // hierarchical grid
      GridPartType& gridPart_; // grid part(view), here the leaf grid the discrete space is build with
 
      const DiscreteFunctionSpaceType& discreteSpace_; // discrete function space
      DiscreteFunctionType residual_;   // initial residual
      DiscreteFunctionType update_;     // distance to solution
 
      DirichletBoundaryFunctionType boundaryValues_;
      DirichletWeightFunctionType boundaryWeight_;
      ConstraintsOperatorType  dirichletConstraints_; // dirichlet boundary constraints
      EmptyConstraintsType     emptyConstraints_;     // do-nothing constraints.
 
      ImplicitOperatorType implicitOperator_; // the affine-linear operator.
      ExplicitOperatorType explicitOperator_;
 
      LinearOperatorType linearOperator_;  // the linear operator (i.e. jacobian of the operator)
 
      const double solverEps_ ; // eps for linear solver
 
      SumModelType sumModel_;
      BulkForcesFunctionType bulkForces_;
      BoundaryFluxFunctionType boundaryFlux_;
 
      mutable int sequence_;            // sequence number
 
      IOTupleType ioTuple_; // tuple with pointers
      DataOutputType *dataOutput_; // data output class
    };
 
    template<class DiscreteFunction,
             class ImplicitModel,
             class ExplicitDataFunction = ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                           typename ImplicitModel::GridPartType>,
             class ExplicitModel = ZeroModel<typename ImplicitModel::FunctionSpaceType,
                                             typename ImplicitModel::GridPartType>,
             class InitialGuessFunction = ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                           typename ImplicitModel::GridPartType> >
    class SplittingFemScheme;
 
    template<class DiscreteFunction,
             class ImplicitModel,
             class ExplicitDataFunction,
             class ExplicitModel,
             class InitialGuessFunction>
    class SplittingFemScheme
      : public SplittingFemSchemeBase<DiscreteFunction, ImplicitModel, ExplicitDataFunction, ExplicitModel, InitialGuessFunction>
    {
      typedef
      SplittingFemSchemeBase<DiscreteFunction, ImplicitModel, ExplicitDataFunction, ExplicitModel, InitialGuessFunction>
      BaseType;
     public:
      typedef DiscreteFunction DiscreteFunctionType;
      typedef ImplicitModel ImplicitModelType;
      typedef ExplicitDataFunction ExplicitDataType;
      typedef ExplicitModel ExplicitModelType;
      typedef InitialGuessFunction InitialGuessType;
 
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const ExplicitDataType& explicitData,
                         const ExplicitModelType& explicitModel,
                         const InitialGuessType& initialGuess,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel, explicitData, explicitModel, initialGuess, name)
      {}
    };
 
    // Instantiate default stuff. This cannot be stuffed in to one
    // class because default arguments must not depend on other
    // arguments.
 
    template<class DiscreteFunction,
             class ImplicitModel,
             class ExplicitDataFunction,
             class ExplicitModel>
    class SplittingFemScheme<DiscreteFunction, ImplicitModel, ExplicitDataFunction, ExplicitModel,
                             ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                              typename ImplicitModel::GridPartType> >
      : public SplittingFemSchemeBase<DiscreteFunction, ImplicitModel, ExplicitDataFunction, ExplicitModel,
                                      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                       typename ImplicitModel::GridPartType> >
 
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
      typedef ImplicitModel ImplicitModelType;
      typedef ExplicitDataFunction ExplicitDataType;
      typedef ExplicitModel ExplicitModelType;
      typedef
      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType, typename ImplicitModel::GridPartType>
      InitialGuessType;
     private:
      typedef
      SplittingFemSchemeBase<DiscreteFunctionType, ImplicitModelType, ExplicitDataType, ExplicitModelType, InitialGuessType>
      BaseType;
     public:
 
      // Normal constructor
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const ExplicitDataType& explicitData,
                         const ExplicitModelType& explicitModel,
                         const InitialGuessType& initialGuess,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel, explicitData, initialGuess, name)
      {}
 
      // Constructor with artificial default argument
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const ExplicitDataType& explicitData,
                         const ExplicitModelType& explicitModel,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel, explicitData, explicitModel,
                   InitialGuessType(solution.space().gridPart()),
                   name)
      {}
    };
 
    template<class DiscreteFunction, class ImplicitModel, class InitialGuess>
    class SplittingFemScheme<DiscreteFunction, ImplicitModel,
                             ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                              typename ImplicitModel::GridPartType>,
                             ZeroModel<typename ImplicitModel::FunctionSpaceType,
                                       typename ImplicitModel::GridPartType>,
                             InitialGuess>
      : public SplittingFemSchemeBase<DiscreteFunction, ImplicitModel,
                                      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                       typename ImplicitModel::GridPartType>,
                                      ZeroModel<typename ImplicitModel::FunctionSpaceType,
                                                typename ImplicitModel::GridPartType>,
                                      InitialGuess>
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
      typedef ImplicitModel ImplicitModelType;
      typedef
      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType, typename ImplicitModel::GridPartType>
      ExplicitDataType;
      typedef
      ZeroModel<typename ImplicitModel::FunctionSpaceType, typename ImplicitModel::GridPartType>
      ExplicitModelType;
      typedef InitialGuess InitialGuessType;
     private:
      typedef
      SplittingFemSchemeBase<DiscreteFunctionType, ImplicitModelType, ExplicitDataType, ExplicitModelType, InitialGuessType>
      BaseType;
     public:
 
      // Normal constructor
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const ExplicitDataType& explicitData,
                         const ExplicitModelType& explicitModel,
                         const InitialGuessType& initialGuess,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel, explicitData, initialGuess, name)
      {}
 
      // Constructor with artificial default arguments
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const InitialGuessType& initialGuess,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel,
                   ExplicitDataType(solution.space().gridPart()),
                   ExplicitModelType(),
                   initialGuess,
                   name)
      {}
    };
 
    template<class DiscreteFunction, class ImplicitModel>
    class SplittingFemScheme<DiscreteFunction, ImplicitModel,
                                      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                       typename ImplicitModel::GridPartType>,
                                      ZeroModel<typename ImplicitModel::FunctionSpaceType,
                                                typename ImplicitModel::GridPartType>,
                                      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                       typename ImplicitModel::GridPartType> >
      : public SplittingFemSchemeBase<DiscreteFunction, ImplicitModel,
                                      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                       typename ImplicitModel::GridPartType>,
                                      ZeroModel<typename ImplicitModel::FunctionSpaceType,
                                                typename ImplicitModel::GridPartType>,
                                      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType,
                                                       typename ImplicitModel::GridPartType> >
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
      typedef ImplicitModel ImplicitModelType;
      typedef
      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType, typename ImplicitModel::GridPartType>
      ExplicitDataType;
      typedef
      ZeroModel<typename ImplicitModel::FunctionSpaceType, typename ImplicitModel::GridPartType>
      ExplicitModelType;
      typedef
      ZeroGridFunction<typename ImplicitModel::FunctionSpaceType, typename ImplicitModel::GridPartType>
      InitialGuessType;
     private:
      typedef
      SplittingFemSchemeBase<DiscreteFunctionType, ImplicitModelType, ExplicitDataType, ExplicitModelType, InitialGuessType>
      BaseType;
     public:
 
      // Normal constructor
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const ExplicitDataType& explicitData,
                         const ExplicitModelType& explicitModel,
                         const InitialGuessType& initialGuess,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel, explicitData, initialGuess, name)
      {}
 
      // Constructor with artificial default arguments
      SplittingFemScheme(DiscreteFunctionType& solution,
                         const ImplicitModelType& implicitModel,
                         const std::string name = "femscheme")
        : BaseType(solution, implicitModel,
                   ExplicitDataType(solution.space().gridPart()),
                   ExplicitModelType(),
                   InitialGuessType(solution.space().gridPart()),
                   name)
      {}
    };
 
    template <class DiscreteFunction, class ImplicitModel, class ExplicitDataFunction, class ExplicitModel>
    void singleSolve(DiscreteFunction& solution,
                     const ImplicitModel& implicitModel,
                     const ExplicitDataFunction& explicitData,
                     const ExplicitModel& explicitModel,
                     const std::string name = "femscheme")
    {
      SplittingFemScheme<DiscreteFunction,ImplicitModel,ExplicitDataFunction,ExplicitModel>
        scheme( solution, implicitModel, explicitData, explicitModel, name );
      scheme.solve(true);
    }
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif //  __DUNE_ACFEM_SPLITTINGFEMSCHEME_HH__