surfmassmatrix.hh

#ifndef SURFACE_MASSMATRIX_HH
#define SURFACE_MASSMATRIX_HH

#if HAVE_DUNE_DISC

/** \file
    \brief Contains a method which assembles the FE mass matrix of a boundary patch of a grid.
*/
#include <dune/grid/common/genericreferenceelements.hh>
#include <dune/grid/common/quadraturerules.hh>

#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>

#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>

#include "boundarypatch.hh"

template <class GridView, int blocksize>
void assembleSurfaceMassMatrix(const BoundaryPatchBase<GridView>& boundary,
                               Dune::BCRSMatrix<Dune::FieldMatrix<double,blocksize,blocksize> >& matrix)
{
    using namespace Dune;

    const int dim = GridView::dimension;
    const GridView& gridView = boundary.gridView();
    const typename GridView::IndexSet& indexSet = gridView.indexSet();

    std::vector<int> globalToLocal;
    boundary.makeGlobalToLocal(globalToLocal);

    // /////////////////////////////////////////////////
    //   Compute matrix occupation structure
    // /////////////////////////////////////////////////

    MatrixIndexSet indices(boundary.numVertices(), boundary.numVertices());

    typename BoundaryPatchBase<GridView>::iterator it    = boundary.begin();
    typename BoundaryPatchBase<GridView>::iterator endIt = boundary.end();

    for (; it!=endIt; ++it) {

        const GenericReferenceElement<double,dim>& refElement = GenericReferenceElements<double, dim>::general(it->inside()->type());
        int nDofs = refElement.size(it->indexInInside(),1,dim);
        
        for (int i=0; i<nDofs; i++) {
            int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
            int globalRow = globalToLocal[indexSet.subIndex(*it->inside(), faceIdxi, dim)];
            
            for (int j=0; j<nDofs; j++) {
                int faceIdxj = refElement.subEntity(it->indexInInside(), 1, j, dim);
                int globalCol = globalToLocal[indexSet.subIndex(*it->inside(), faceIdxj, dim)];
                
                indices.add(globalRow,globalCol);
                
            }
            
        } 

    }


    // /////////////////////////////////////////////////
    //   Compute matrix values
    // /////////////////////////////////////////////////

    indices.exportIdx(matrix);
    matrix = 0;

    for (it=boundary.begin(); it!=endIt; ++it) {

        const typename BoundaryPatchBase<GridView>::iterator::Intersection::Geometry& segmentGeometry = it->geometry();

            const GenericReferenceElement<double,dim>& refElement = GenericReferenceElements<double, dim>::general(it->inside()->type());
            int nDofs = refElement.size(it->indexInInside(),1,dim);

            // Local matrix
            Matrix<FieldMatrix<double,1,1> > scalarMat(nDofs, nDofs);
            scalarMat = 0;

            int quadOrder = 2*dim;

            // Get quadrature rule
            const QuadratureRule<double, dim-1>& quad = QuadratureRules<double, dim-1>::rule(segmentGeometry.type(), quadOrder);

            // Get set of shape functions on this segment
            // TODO: We should really use the new shape functions here, since we use the new
            // local indices below!!! This only works since intervall, triangle and square nodes numbering
            // remains the same.
            const typename LagrangeShapeFunctionSetContainer<double,double,dim-1>::value_type& sfs
                = LagrangeShapeFunctions<double,double,dim-1>::general(segmentGeometry.type(),1);

            /* Loop over all integration points */
            for (int ip=0; ip<quad.size(); ip++) {
                
                // Local position of the quadrature point
                const FieldVector<double,dim-1>& quadPos = quad[ip].position();
                
                const double integrationElement = segmentGeometry.integrationElement(quadPos);

                // Evaluate base functions
                double v[nDofs];

                for(int i=0; i<nDofs; i++)
                    v[i] = sfs[i].evaluateFunction(0,quadPos);

                for(int i=0; i<nDofs; i++) 
                    for (int j=0; j<=i; j++ ) 
                        scalarMat[i][j] += ( v[i] * v[j] ) * quad[ip].weight() * integrationElement;
            
            }

            // Make the matrix symmetric
            for(int i=0; i<nDofs; i++)
                for (int j=nDofs-1; j>i; j--)
                    scalarMat[i][j] = scalarMat[j][i];

            for (int i=0; i<nDofs; i++) {
                int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
                int globalRow = globalToLocal[indexSet.subIndex(*it->inside(), faceIdxi, dim)];
                
                for (int j=0; j<nDofs; j++) {
                    int faceIdxj = refElement.subEntity(it->indexInInside(), 1, j, dim);
                    int globalCol = globalToLocal[indexSet.subIndex(*it->inside(), faceIdxj, dim)];

                    for (int k=0; k<blocksize; k++)
                        matrix[globalRow][globalCol][k][k] += scalarMat[i][j];

                }

            } 

    }

}

#else 
    #warning surfmassmatrix.hh was included but dune-disc is not present !
#endif // end of #if HAVE_DUNE_DISC

#endif