functionalassembler.hh

#ifndef FUNCTIONAL_ASSEMBLER_HH
#define FUNCTIONAL_ASSEMBLER_HH

#include <dune/istl/bvector.hh>

#include "dune/fufem/functionspacebases/functionspacebasis.hh"

//! Generic global assembler for functionals on a gridview
template <class TrialBasis>
class FunctionalAssembler
{
    private:
        typedef typename TrialBasis::GridView GridView;

    public:
        //! create assembler for gridview
        FunctionalAssembler(const TrialBasis& tBasis) :
            tBasis_(tBasis)
        {}

        /** \brief Assemble
         *
         * \param localAssembler local assembler
         * \param[out] b target vector
         * \param initializeVector If this is set the output vector is
         * set to the correct length and initialized to zero before
         * assembly. Otherwise the assembled values are just added to
         * the vector.
         */
        template <class LocalFunctionalAssemblerType, class GlobalVectorType>
        void assemble(LocalFunctionalAssemblerType& localAssembler, GlobalVectorType& b, bool initializeVector=true) const
        {
            typedef typename GridView::template Codim<0>::Iterator ElementIterator;
            typedef typename LocalFunctionalAssemblerType::LocalVector LocalVector;
            typedef typename TrialBasis::LinearCombination LinearCombination;

            int rows = tBasis_.size();

            if (initializeVector)
            {
                b.resize(rows);
                b=0.0;
            }

            ElementIterator it = tBasis_.getGridView().template begin<0>();
            ElementIterator end = tBasis_.getGridView().template end<0>();
            for (; it != end; ++it)
            {
                // get shape functions
                const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(*it);

                LocalVector localB(tFE.localBasis().size());
                localAssembler.assemble(*it, localB, tFE);

                for (size_t i=0; i<tFE.localBasis().size(); ++i)
                {
                    int idx = tBasis_.index(*it, i);
                    const LinearCombination& constraints = tBasis_.constraints(idx);
                    bool isConstrained = tBasis_.isConstrained(idx);
                    if (isConstrained)
                    {
                        for (size_t w=0; w<constraints.size(); ++w)
                            b[constraints[w].index].axpy(constraints[w].factor, localB[i]);
                    }
                    else
                        b[idx] += localB[i];
                }
            }
            return;
        }

    private:
        const TrialBasis& tBasis_;
};

#endif