Skip to content
Snippets Groups Projects
Forked from agnumpde / dune-fufem
1862 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
gradientassembler.hh 3.39 KiB
#ifndef GRADIENT_ASSEMBLER_HH
#define GRADIENT_ASSEMBLER_HH


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

#include <dune/grid/common/quadraturerules.hh>
#include <dune/grid/common/genericreferenceelements.hh>

#include <dune/istl/matrix.hh>

#include <dune/ag-common/functions/localbasisderivativefunction.hh>
#include <dune/ag-common/assemblers/localassembler.hh>


//! \todo Please doc me !
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
class GradientAssembler
    : public LinearLocalAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix<double,GridType::dimension,1> >
{

        static const int dim = GridType::dimension;

    public:
        typedef typename Dune::FieldMatrix<double,dim,1> T;

        typedef typename GridType::template Codim<0>::Entity::EntityPointer ElementPointer;
        typedef typename LocalAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >::BoolMatrix BoolMatrix;
        typedef typename LocalAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >::LocalMatrix LocalMatrix;

        void indices(const ElementPointer& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
        {
            isNonZero = true;
            return;
        }

        void assemble(const ElementPointer& element, LocalMatrix& localMatrix
                  , const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
        {
            typedef typename Dune::template FieldVector<double,dim> FVdim;
            typedef typename Dune::template FieldMatrix<double,dim,dim> FMdimdim;
            typedef typename AnsatzLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
            typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::DomainType DomainType;
            typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::RangeType RangeType;

            localMatrix = 0.0;

            // compute the element barycenter in local coordinates
            const FVdim pos = Dune::GenericReferenceElements<double,dim>::general(element->type()).position(0,0);

            // get transposed inverse of Jacobian of transformation
            const FMdimdim& invJacobian = element->geometry().jacobianInverseTransposed(pos);

            typedef typename Dune::LocalFiniteElementFunctionBase<TrialLocalFE>::type FunctionBaseClass;

            LocalBasisDerivativeFunction<AnsatzLocalFE, FunctionBaseClass> derivative(aFE.localBasis());

            std::vector<JacobianType> referenceGradients(tFE.localBasis().size());

            std::vector<typename RangeType::field_type> partialDerivatives(tFE.localBasis().size());

            for (int i = 0; i < aFE.localBasis().size(); ++i)
            {
                // interpolate all partial derivatives of ansatz function by test functions
                derivative.setIndex(i);
                for (int j = 0; j < dim; ++j)
                {
                    derivative.setComponent(j);
                    tFE.localInterpolation().interpolate(derivative, partialDerivatives);

                    for (int k = 0; k < tFE.localBasis().size(); ++k)
                        referenceGradients[k][0][j] = partialDerivatives[k];
                }

                // transform gradients
                for (int k = 0; k < tFE.localBasis().size(); ++k)
                    invJacobian.mv(referenceGradients[k][0], localMatrix[k][i]);
            }
        }

};


#endif