boundarymassassembler.hh

#ifndef BOUNDARY_MASS_ASSEMBLER_HH
#define BOUNDARY_MASS_ASSEMBLER_HH


#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <dune/istl/matrix.hh>

#include "dune/fufem/staticmatrixtools.hh"
#include "dune/fufem/quadraturerules/quadraturerulecache.hh"

#include "dune/fufem/assemblers/localoperatorassembler.hh"
#include "dune/fufem/referenceelementhelper.hh"

//** \brief Local mass assembler **//
template <class GridType, class BoundaryPatchType, class TrialLocalFE, class AnsatzLocalFE, class T=Dune::FieldMatrix<double,1,1> >
class BoundaryMassAssembler : public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >
{
    private:
        static const int dim = GridType::dimension;

    public:
        typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE ,T >::Element Element;
        typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE ,T >::BoolMatrix BoolMatrix;
        typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE ,T >::LocalMatrix LocalMatrix;
        typedef typename BoundaryPatchType::iterator::Intersection Intersection;

        BoundaryMassAssembler(const BoundaryPatchType& boundaryPatch, int quadOrder=2):
            boundaryPatch_(boundaryPatch),
            quadOrder_(quadOrder)
        {}

        void indices(const Element& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
        {
            if (boundaryPatch_.containsFaceOf(element))
                isNonZero = true;
            else
                isNonZero = false;
        }

        void assemble(const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
        {
            typedef typename BoundaryPatchType::GridView BoundaryPatchGridView;
            typedef typename BoundaryPatchType::GridView::IntersectionIterator IntersectionIterator;

            // Make sure we got suitable shape functions
            assert(tFE.type() == element.type());
            assert(aFE.type() == element.type());

            localMatrix = 0.0;

            if (not(boundaryPatch_.containsFaceOf(element)))
                return;

            const BoundaryPatchGridView& gridView = boundaryPatch_.gridView();

            IntersectionIterator it = gridView.ibegin(element);
            IntersectionIterator end = gridView.iend(element);
            for(; it != end; ++it)
            {
                const Intersection& i = *it;
                if (boundaryPatch_.contains(i))
                    assembleOnIntersection(i, element, localMatrix, tFE, aFE);
            }
        }

        void assembleOnIntersection(const Intersection& intersection, const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
        {
            typedef typename Dune::template FieldVector<double,dim> ElementLocalCoordinate;
            typedef typename Dune::template FieldVector<double,dim-1> IntersectionLocalCoordinate;
            typedef typename Intersection::Geometry IntersectionGeometry;
            typedef typename Intersection::LocalGeometry IntersectionGeometryInElement;
            typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::RangeType RangeType;

            // Make sure we got suitable shape functions
            assert(tFE.type() == element.type());
            assert(aFE.type() == element.type());

            const int tFEsize = tFE.localBasis().size();
            const int aFEsize = aFE.localBasis().size();

            localMatrix = 0.0;

            // get quadrature rule
//            const Dune::template QuadratureRule<double, dim>& quad = Dune::template QuadratureRules<double, dim>::rule(element.type(), quadOrder_);
            const Dune::template QuadratureRule<double, dim-1>& quad = QuadratureRuleCache<double, dim-1>::rule(intersection.type(), quadOrder_, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );

            const IntersectionGeometry intersectionGeometry = intersection.geometry();
            const IntersectionGeometryInElement intersectionGeometryInElement = intersection.geometryInInside();

            // store values of shape functions
            std::vector<RangeType> tFEvalues(tFEsize);
            std::vector<RangeType> aFEvalues(aFEsize);

            // determine if local dof is located on intersection
            typename std::vector<char> intersectionContainsTFE(tFEsize);
            for (int i = 0; i < tFEsize; ++i)
                intersectionContainsTFE[i] = intersectionContainsLocalKey(element.type(), intersection, tFE.localCoefficients().localKey(i));
            typename std::vector<char> intersectionContainsAFE(aFEsize);
            for (int i = 0; i < aFEsize; ++i)
                intersectionContainsAFE[i] = intersectionContainsLocalKey(element.type(), intersection, aFE.localCoefficients().localKey(i));

            // loop over quadrature points
            for (size_t pt=0; pt < quad.size(); ++pt)
            {
                // get quadrature point
                const IntersectionLocalCoordinate& intersectionQuadPos = quad[pt].position();
                ElementLocalCoordinate quadPos = intersectionGeometryInElement.global(intersectionQuadPos);

                // get integration factor
                const double integrationElement = intersectionGeometry.integrationElement(intersectionQuadPos);

                // evaluate basis functions
                tFE.localBasis().evaluateFunction(quadPos, tFEvalues);
                aFE.localBasis().evaluateFunction(quadPos, aFEvalues);

                // compute matrix entries
                double z = quad[pt].weight() * integrationElement;
                for(int i=0; i<tFEsize; ++i)
                {
                    if (intersectionContainsTFE[i])
                    {
                        double zi = tFEvalues[i]*z;

                        for (int j=0; j<aFEsize; ++j)
                        {
                            if (intersectionContainsAFE[j])
                                StaticMatrix::addToDiagonal(localMatrix[i][j], aFEvalues[j] * zi);
                        }
                    }
                }
            }
        }
    protected:

        static bool intersectionContainsLocalKey(const Dune::GeometryType& gt, const Intersection& intersection, const typename Dune::LocalKey& localKey)
        {
            return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(gt, intersection.indexInInside(), 1, localKey.subEntity(), localKey.codim());
        }


        const BoundaryPatchType& boundaryPatch_;

        const int quadOrder_;
};


#endif