l2projection.hh

#ifndef __DUNE_ACFEM_L2PROJECTION_HH__
#define __DUNE_ACFEM_L2PROJECTION_HH__
 
#include "functionals/linearfunctional.hh"
#include "../models/basicmodels.hh"
#include "../models/modelexpression.hh"
#include "../algorithms/ellipticfemscheme.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunction, class Functional>
    class L2ProjectionOperator
    {
     public:
      typedef Functional FunctionalType;
      typedef typename FunctionalType::TraitsType TraitsType;
      typedef typename FunctionalType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef DiscreteFunction DiscreteFunctionType;
      typedef
      DiscreteLinearFunctional<DiscreteFunctionSpaceType, TraitsType>
      FunctionalInterfaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
     public:
 
      void operator()(const FunctionalInterfaceType& phi,
                      DiscreteFunctionType& result)
      {
        // skip the initialization, just do a quick loop over "result"
        // in order to determine that all values are finite
        auto dend(result.dend());
        for (auto it = result.dbegin(); it != dend; ++it) {
          if (!std::isfinite(*it)) {
            result.clear();
            break;
          }
        }
 
        // In principle it should be fine to stuff a MassModel into the
        // EllipticFemScheme and only call the method scheme.solve() ...
        //
        // Let's try. We do not maintain any state, the constructor is
        // trivial. After all, as all default implementations in
        // principle should reduce to nothing, this should be fairly
        // efficient.
 
        auto model = massModel(result) - phi;
        typedef
          EllipticFemScheme<DiscreteFunctionType, decltype(model)>
          SolverSchemeType;
 
        SolverSchemeType solverScheme(result, model);
 
        solverScheme.solve();
 
        // And that should be it ...
      }
    };
 
    template<class FunctionSpace, class GridFunction, class DiscreteFunction>
    static inline
    void L2Projection(const Fem::Function<FunctionSpace, GridFunction>& fct,
                      DiscreteFunction& result)
    {
      typedef DiscreteFunction DiscreteFunctionType;
      typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef L2InnerProductFunctional<DiscreteFunctionSpaceType, GridFunction> FunctionalType;
      typedef L2ProjectionOperator<DiscreteFunctionType, FunctionalType> L2ProjectionType;
 
      FunctionalType functional(result.space(), fct);
      L2ProjectionType projection;
 
      projection(functional, result);
    }
 
    template<class Traits, class DiscreteFunction>
    static inline
    void L2Projection(const DiscreteLinearFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType, Traits>& phi,
                      DiscreteFunction& result)
    {
      typedef typename Traits::GlobalFunctionalType FunctionalType;
      typedef DiscreteFunction DiscreteFunctionType;
      typedef L2ProjectionOperator<DiscreteFunctionType, FunctionalType> L2ProjectionType;
 
      const FunctionalType& functional(static_cast<const FunctionalType&>(phi));
      L2ProjectionType projection;
 
      projection(functional, result);
    }
 
 
 
  } // ACFem::
 
} // Dune::
 
#endif // __DUNE_ACFEM_L2PROJECTION_HH__