Skip to content
Snippets Groups Projects
Select Git revision
  • adb5aff6f0af8167c6a7592f7c75ec0f17b4524d
  • master default protected
  • releases/2.6-1
  • releases/2.5-1
  • reduce_contact_patches
  • ansgar-2.4-compat
  • feature/large_deformation_contact
  • releases/2.4 protected
  • releases/2.3-2 protected
  • releases/2.3-1 protected
  • releases/2.2-1 protected
  • releases/2.1-1 protected
  • releases/2.0-1 protected
  • subversion->git
14 results

gapfunction.cc

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    gapfunction.cc 36.90 KiB
    // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
    // vi: set ts=8 sw=4 et sts=4:
    #include <dune/fufem/arithmetic.hh>
    #include <dune/fufem/staticmatrixtools.hh>
    #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
    #include <dune/fufem/boundarywriter.hh>
    #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
    #include <dune/fufem/functions/basisgridfunction.hh>
    
    #include <dune/localfunctions/lagrange/p1.hh>
    #include <dune/localfunctions/dualmortarbasis/dualpq1factory.hh>
    #include <dune/grid-glue/adapter/rangegenerators.hh>
    #include <dune/contact/assemblers/linearisedaveragednormalsassembler.hh>
    #include <dune/contact/assemblers/linearisedintersectionassembler.hh>
    
    namespace Dune {
    namespace Contact {
    
    template <class GridType, class field_type>
    void GapFunction<GridType,field_type>::computeMergedGrid()
    {
        // compute contact mapping using GridGlue
        const auto& patch0 = *patches_[0];
        const auto& gridView0 = patch0.gridView();
        const auto& gridView1 = patches_[1]->gridView();
    
        using Element = typename LeafGridView::template Codim<0>::Entity;
    
        auto desc0 = [&] (const Element& e, unsigned int face) {
            return patches_[0]->contains(e,face);
        };
    
        auto desc1 = [&] (const Element& e, unsigned int face) {
            return patches_[1]->contains(e,face);
        };
    
    
        glue_ = std::make_shared<GlueType>(std::make_shared<Extractor>(gridView0,desc0),
                                           std::make_shared<Extractor>(gridView1,desc1),
                                           merger_);
        glue_->build();
    
        reducedNonmortar_.setup(gridView0);
        reducedMortar_.setup(gridView1);
    
        // Restrict the hasObstacle fields to the part of the nonmortar boundary
        // where we could actually create a contact mapping
    
        // Get all fine grid boundary segments that are totally covered by the grid-glue segments
        using Tuple = std::pair<int,int>;
        std::map<Tuple, field_type> coveredArea, fullArea;
    
        const auto& indexSet = gridView0.indexSet();
    
        // initialize with area of boundary faces
        for (const auto& bIt : patch0)
            coveredArea[Tuple(indexSet.index(bIt.inside()),bIt.indexInInside())] = bIt.geometry().volume();
    
        // save areas
        fullArea = coveredArea;
    
        // subtract the remote intersection areas to find out which are totally covered
        for (const auto& rIs : intersections(*glue_))
            coveredArea[Tuple(indexSet.index(rIs.inside()),rIs.indexInInside())] -=
                    rIs.geometry().volume();
    
        // add all fine grid faces that are totally covered by the contact mapping
        for (const auto& bIt : patch0) {
            const auto inside = bIt.inside();
            if(coveredArea[Tuple(indexSet.index(inside),bIt.indexInInside())]/
                    fullArea[Tuple(indexSet.index(inside),bIt.indexInInside())] < 0.001)
                reducedNonmortar_.addFace(inside, bIt.indexInInside());
        }
    
    
        for (const auto& rIs : intersections(*glue_))
            if (reducedNonmortar_.contains(rIs.inside(),rIs.indexInInside()))
                reducedMortar_.addFace(rIs.outside(),rIs.indexInOutside());
    
        std::cout<<"Reduced Nonmortar contains "<<reducedNonmortar_.numFaces()<<" of "
                <<patch0.numFaces()<<" original nonmortar faces "<<std::endl;
    }
    
    template <class GridType, class field_type>
    Dune::BlockVector<Dune::FieldVector<field_type,1> >
    GapFunction<GridType,field_type>::operator()(const VectorType& nonmortarDisp,
                                                         const VectorType& mortarDisp, bool glue)
    {
        if (glue)
            computeMergedGrid(nonmortarDisp,mortarDisp);
    
        // get index sets
        const auto& gridView0 = complex_.grid("0").leafGridView();
        const auto& indexSet0 = gridView0.indexSet();
    
        // cache for the dual functions on the boundary
        typedef Dune::DualPQ1LocalFiniteElementCache<ctype, field_type, dim,true> DualFECache;
        DualFECache dualCache;
    
        std::vector<int> globalToLocal;
        reducedNonmortar_.makeGlobalToLocal(globalToLocal);
    
        // assemble weak obstacles
        RangeType weakObstacle(reducedNonmortar_.numVertices());
        weakObstacle = 0;
    
        Dune::GeometryType remoteGt = Dune::GeometryTypes::simplex(dim-1);
    
        // loop over all intersections
        std::vector<CoordinateType> avNormals = patches_[1]->getNormals();
        using Basis = P1NodalBasis<typename DeformedGridType::LeafGridView>;
        Basis p1Basis(complex_.grid("1").leafGridView());
        BasisGridFunction<Basis, std::vector<CoordinateType> >  mortarNormal(p1Basis, avNormals);
    
        for (const auto& rIs : intersections(*glue_)) {
    
            // types of the elements supporting the boundary segments in question
            const auto& inside = rIs.inside();
    
            if (!reducedNonmortar_.contains(inside,rIs.indexInInside()))
                continue;
    
            const auto& outside = rIs.outside();
    
            Dune::GeometryType nmElemType = inside.type();
    
            const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nmElemType);
            const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
    
            Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
            Dune::GeometryType mFaceType  = targetRefElement.type(rIs.indexInOutside(),1);
    
            // Select a quadrature rule
            // 2 in 2d and for integration over triangles in 3d.  If one (or both) of the two faces involved
            // are quadrilaterals, then the quad order has to be risen to 3 (4).
            /****************************************** CHANGE BACK TO 2 ******************************/
            int quadOrder = 4 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
            const auto& quad = QuadratureRuleCache<ctype, dim-1>::rule(remoteGt, quadOrder, false);
    
            const auto& dualFiniteElement = dualCache.get(nmElemType);
            using LocalBasisType = typename DualFECache::FiniteElementType::Traits::LocalBasisType;
    
            std::vector<typename LocalBasisType::Traits::RangeType> dualQuadValues;
    
            const auto& rGeom = rIs.geometry();
            const auto& rGeomOutside = rIs.geometryOutside();
            const auto& rGeomInInside = rIs.geometryInInside();
            const auto& rGeomInOutside = rIs.geometryInOutside();
    
            std::vector<int> nmFaceIdx(domainRefElement.size(rIs.indexInInside(),1,dim));
            for (size_t i=0; i<nmFaceIdx.size(); i++)
              nmFaceIdx[i] = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
    
            for (size_t l=0; l<quad.size(); l++) {
    
                // compute integration element of overlap
                ctype integrationElement = rGeom.integrationElement(quad[l].position());
    
                // quadrature point positions on the reference element
                CoordinateType nonmortarQuadPos = rGeomInInside.global(quad[l].position());
    
                // The current quadrature point in world coordinates
                CoordinateType nonmortarQpWorld = rGeom.global(quad[l].position());
                CoordinateType mortarQpWorld    = rGeomOutside.global(quad[l].position());;
    
                CoordinateType normalQp;
                mortarNormal.evaluateLocal(outside, rGeomInOutside.global(quad[l].position()), normalQp);
                normalQp /= normalQp.two_norm();
                auto gap = (nonmortarQpWorld  - mortarQpWorld)*normalQp;
    
                //evaluate all dual functions at the quadrature point
                dualFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);
    
                // loop over all dual shape functions
                for (size_t j=0; j<nmFaceIdx.size(); j++) {
    
                    int globalDomainIdx = indexSet0.subIndex(inside,nmFaceIdx[j],dim);
                    int rowIdx = globalToLocal[globalDomainIdx];
    
                    // if the vertex is not contained in the restricted contact patch then dismiss it
                    weakObstacle[rowIdx][0] += integrationElement * quad[l].weight()* dualQuadValues[nmFaceIdx[j]] * gap;
                }
            }
        }
        return weakObstacle;
    }
    
    template <class GridType, class field_type>
    void GapFunction<GridType,field_type>::firstDerivative(const VectorType& nonmortarDisp,
            const VectorType& mortarDisp,
            JacobianType& linearMat, bool glue)
    {
        std::cout<<"Compute exact linearization of the mortar constraints!\n";
        namespace MV = Dune::MatrixVector;
    
        // check if merged grids need to be computed
        if (glue)
            computeMergedGrid(nonmortarDisp,mortarDisp);
    
        const auto& indexSet0 = complex_.grid("0").leafIndexSet();
        const auto& indexSet1 = complex_.grid("1").leafIndexSet();
        int size0 = indexSet0.size(dim);
        int size1 = indexSet1.size(dim);
    
        // cache for the dual functions on the element
        typedef Dune::DualPQ1LocalFiniteElementCache<ctype, field_type, dim,true> DualFECache;
        DualFECache dualCache;
    
        // cache for the Lagrange shapefunctions element
        typedef Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> LagrangeFECache;
        LagrangeFECache femCache;
    
        using DualFaceFECache = DualPQ1LocalFiniteElementCache<ctype, field_type, dim-1>;
        DualFaceFECache dualFaceCache;
    
        // cache for Lagrange shape functions on the faces we need them to compute
        // the linearisation of the local geometries of the remote intersecitons more stable
        using LagrangeFaceFECache = PQkLocalFiniteElementCache<ctype, field_type, dim-1, 1>;
        LagrangeFaceFECache femFaceCache;
    
        // Create mapping from the global set of block dofs to the ones on the contact boundary
        std::vector<int> globalToLocal;
        reducedNonmortar_.makeGlobalToLocal(globalToLocal);
    
        std::vector<int> globalToLocalM;
        reducedMortar_.makeGlobalToLocal(globalToLocalM);
    
        int nNMvert = reducedNonmortar_.numVertices();
    
        // get the nodal averaged unit normals
        std::vector<CoordinateType> avNormalsM = patches_[1]->getNormals();
    
        using Basis = P1NodalBasis<typename DeformedGridType::LeafGridView>;
        Basis p1Basis(complex_.grid("1").leafGridView());
        BasisGridFunction<Basis, std::vector<CoordinateType> >  mortarNormal(p1Basis, avNormalsM);
    
        // compute derivative of the averaged unit normals w.r.t the mortar coefficients
        MatrixType avNormMortarDeriv;
        LinearisedAveragedNormalsAssembler::assemble(reducedMortar_, *patches_[1], avNormMortarDeriv);
    
        // compute derivative of the remote vertices w.r.t the nonmortar coefficients
        MatrixType remoteVerticesDeriv;
        LinearisedIntersectionAssembler<GlueType> remoteAssembler(glue_);
        remoteAssembler.assemble(reducedNonmortar_, *patches_[0],
                                 reducedMortar_, *patches_[1],
                                 merger_->getOverlap(), std::cos(merger_->minNormalAngle()),
                                 remoteVerticesDeriv, avNormMortarDeriv);
    
        // ///////////////////////////////////////////////////////////
        //   Get the occupation structure for the linearization matrix
        // ///////////////////////////////////////////////////////////
    
        const auto& remoteIntersectionIndexSet = glue_->indexSet();
    
        Dune::MatrixIndexSet mortarIndices(nNMvert, size0+size1);
    
        // loop over all intersections
        for (const auto& rIs : intersections(*glue_)) {
    
            const auto& inside = rIs.inside();
            const auto& outside = rIs.outside();
    
            // if intersection is not contained then skip it
            if (!reducedNonmortar_.contains(inside,rIs.indexInInside()))
                continue;
    
            const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
            const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
    
            int nNonMortVert = domainRefElement.size(dim);
            int nMortVert = targetRefElement.size(dim);
    
            for (int j=0; j<nNonMortVert; j++) {
    
                int globalDomainIdx = indexSet0.subIndex(inside,j,dim);
                int rowIdx = globalToLocal[globalDomainIdx];
    
                // if the vertex is not contained in the restricted contact boundary then dismiss it
                if (rowIdx == -1)
                    continue;
    
                mortarIndices.add(rowIdx,globalDomainIdx);
    
                // coupling with nonmortar dofs
                for (int k=j+1; k<nNonMortVert; k++) {
    
                    int colIdxk = indexSet0.subIndex(inside,k,dim);
                    int rowIdxk = globalToLocal[colIdxk];
    
                    // if the vertex is not contained in the restricted contact boundary then dismiss it
                    if (rowIdxk != -1) {
                        mortarIndices.add(rowIdx, colIdxk);
                        mortarIndices.add(rowIdxk, globalDomainIdx);
                    }
                }
    
                // coupling with mortar dofs
                for (int k=0; k<nMortVert; k++) {
    
                    int colIdx = indexSet1.subIndex(outside,k,dim);
                    // if the vertex is not contained in the restricted contact boundary then dismiss it
                    if (!reducedMortar_.containsVertex(colIdx))
                        continue;
    
                    mortarIndices.add(rowIdx, size0 + colIdx);
    
                    // add all indices of the linearised mortar normals
                    const auto& cEndIt = avNormMortarDeriv[globalToLocalM[colIdx]].end();
                    for(auto cIt = avNormMortarDeriv[globalToLocalM[colIdx]].begin(); cEndIt != cIt; cIt++)
                      mortarIndices.add(rowIdx,size0 + cIt.index());
                }
    
                // add all indices of the linearised remote vertices
                int rIdx = remoteIntersectionIndexSet.index(rIs)*dim;
                for (int k=0; k<dim; k++) {
                    const auto& cEndIt = remoteVerticesDeriv[rIdx+k].end();
                    for(auto cIt = remoteVerticesDeriv[rIdx+k].begin(); cEndIt != cIt; cIt++)
                        mortarIndices.add(rowIdx,cIt.index());
                }
            }
        }
        mortarIndices.exportIdx(linearMat);
        linearMat = 0;
    
        // the remote geometry type is always a simplex
        Dune::GeometryType remoteGt = Dune::GeometryTypes::simplex(dim-1);
    
        // get local basis for the remote intersection type, this is always a simplex
        // cache for the remote shapefunctions on the boundary
        Dune::P1LocalBasis<ctype,field_type,dim-1> remoteLfe;
        std::vector<Dune::FieldVector<field_type,1> > remoteQuadValues;
    
        // loop over remote intersections
        for (const auto& rIs : intersections(*glue_)) {
    
            const auto& inside = rIs.inside();
    
            // if intersection is not contained then skip it
            if (!reducedNonmortar_.contains(inside,rIs.indexInInside()))
                continue;
    
            const auto& outside = rIs.outside();
    
            // TODO Adjust quadrature for non-simplices
            // 4 is required for simplices
            int quadOrder = 4;
            const auto& quad = QuadratureRuleCache<ctype, dim-1>::rule(remoteGt, quadOrder,false);
    
            // get local basis functions for the dual and mortar elements
            const auto& nonmortarLfe = femCache.get(inside.type()).localBasis();
            const auto& dualLfe = dualCache.get(inside.type()).localBasis();
            const auto& mortarLfe = femCache.get(outside.type()).localBasis();
    
            std::vector<Dune::FieldMatrix<field_type,1,dim> > nmRefJacobians, mRefJacobians, dualQuadRefJacobians;
            std::vector<Dune::FieldVector<field_type,1> > mQuadValues, nmQuadValues, dualQuadValues;
    
            // create a mapping to check which element basis functions correspond to the face
            const auto& ref0 = Dune::ReferenceElements<ctype,dim>::general(inside.type());
            std::vector<int> nmIdx(dim),nmglobalIdx(dim);
            for (int i=0; i<ref0.size(rIs.indexInInside(),1,dim); i++) {
                    nmIdx[i]=ref0.subEntity(rIs.indexInInside(),1,i,dim);
                    nmglobalIdx[i]=indexSet0.subIndex(inside,nmIdx[i],dim);
            }
    
            const auto& ref1 = Dune::ReferenceElements<ctype,dim>::general(outside.type());
            std::vector<int> mIdx(dim),mglobalIdx(dim);
            for (int i=0; i<ref1.size(rIs.indexInOutside(),1,dim); i++) {
                    mIdx[i]=ref1.subEntity(rIs.indexInOutside(),1,i,dim);
                    mglobalIdx[i]=size0+indexSet1.subIndex(outside,mIdx[i],dim);
            }
    
            // get finite element of the faces
            const auto& nmFaceLfe = femFaceCache.get(ref0.type(rIs.indexInInside(),1)).localBasis();
            const auto& dualFaceLfe = dualFaceCache.get(ref0.type(rIs.indexInInside(),1)).localBasis();
            const auto& mFaceLfe = femFaceCache.get(ref1.type(rIs.indexInOutside(),1)).localBasis();
    
            // compute local geometries for the faces
            typename DeformedGridType::LeafGridView::Intersection nmIntersection;
            for (auto&& is0 : intersections(complex_.grid("0").leafGridView(), inside))
              if (is0.indexInInside()==rIs.indexInInside()) {
                nmIntersection = std::move(is0);
                break;
              }
    
            typename DeformedGridType::LeafGridView::Intersection mIntersection;
            for (auto&& is1 : intersections(complex_.grid("1").leafGridView(), outside))
              if (is1.indexInInside()==rIs.indexInOutside()) {
                mIntersection = std::move(is1);
                break;
              }
    
            const auto& nmLocalFaceGeometry = nmIntersection.geometryInInside();
            const auto& mLocalFaceGeometry = mIntersection.geometryInInside();
    
            // global geometry of the inside and outside element
            const auto& nmElemGeom = inside.geometry();
            const auto& mElemGeom = outside.geometry();
            // geometries of the remote intersection
            const auto& rGlobGeom = rIs.geometry();
            const auto& rGeomOut = rIs.geometryInOutside();
            const auto& rGeomIn = rIs.geometryInInside();
    
            CoordinateVectorType nmCorners(dim), mCorners(dim);
            for (size_t i=0; i<dim; i++) {
                nmCorners[i] = nmElemGeom.corner(nmIdx[i]);
                mCorners[i] = mElemGeom.corner(mIdx[i]);
            }
    
            // By now simplices are hardwired,
            // still use these names to keep in mind what we are looping over
            size_t nmSize = nmCorners.size();
            size_t mSize = mCorners.size();
    
            // compute mortar triangle
            std::vector<CoordinateType> nmEdges(dim);
            for (size_t i=0; i<dim; i++)
                nmEdges[i] = nmCorners[(i+2)%dim]-nmCorners[(i+1)%dim];
            CoordinateType nmNormal;
            if (dim==3) {
              // we want 1-0 cross 2-0
              nmNormal = MV::crossProduct(nmEdges[1],nmEdges[2]);
            } else if (dim==2)
              nmNormal = {{nmEdges[1][1], nmEdges[0][0]}};
    
            // compute linearisation of the projected nonmortar vertices
            std::array<MatrixType,dim> linVertex;
            CoordinateVectorType projVertex(dim);
    
            for (size_t i=0; i<mSize;i++)
                AnalyticLocalIntersectionAssembler<GlueType,MatrixType>::linearisedForwardProjection(
                      linVertex[i],projVertex[i],size0, 0, size0+size1, mglobalIdx[i]-size0, mCorners[i],
                      avNormalsM[mglobalIdx[i]-size0],avNormMortarDeriv[globalToLocalM[mglobalIdx[i]-size0]],
                    nmglobalIdx, nmCorners[0], nmEdges,nmNormal);
    
            ////////////////////////////////////////////////////////////////////////
            //  Derivation of the integration element w.r.t nonmortar coefficients
            ///////////////////////////////////////////////////////////////////////
    
            int realRemoteIdx = remoteIntersectionIndexSet.index(rIs);
            int rIdx = realRemoteIdx*dim;
    
            // the integration element is assumed to be 2*remoteVolume
            // it is given by the norm of the cross-product of the edges
            // thus it involves the linearisation of the remote corners
            Dune::MatrixIndexSet areaIndexSet(1,size0+size1);
            for (int i=0; i<dim; i++) {
                const auto& row = remoteVerticesDeriv[rIdx+i];
                const auto& cEndIt = row.end();
                for (auto cIt=row.begin(); cIt != cEndIt; cIt++)
                    areaIndexSet.add(0,cIt.index());
            }
            JacobianType linIntElem;
            areaIndexSet.exportIdx(linIntElem);
            linIntElem=0;
    
            // make vector containing the edge vectors
            VectorType edge(dim-1);
            for (size_t i=0; i<edge.size(); i++)
                edge[i] = rGlobGeom.corner(i+1) - rGlobGeom.corner(0);
    
            // a normal of the remote intersection
            CoordinateType n;
            if (dim ==3)
               n = MV::crossProduct(edge[0],edge[1]);
            else if (dim==2)
               n = edge[0]; // this is not the normal
            n /= n.two_norm();
    
            // terms appear in the linearisation
            CoordinateVectorType nTimesEdge(dim);
            if (dim==3) {
                nTimesEdge[2] = MV::crossProduct(n,edge[0]);
                nTimesEdge[1] = MV::crossProduct(edge[1],n);
                nTimesEdge[0] = 0; nTimesEdge[0] -=nTimesEdge[2]+nTimesEdge[1];
            } else if (dim==2) {
                nTimesEdge[0] = 0; nTimesEdge[0].axpy(-1,n);
                nTimesEdge[1] = 0; nTimesEdge[1].axpy(1,n);
            }
    
            for (int i=0; i<dim; i++) {
    
                const auto& cEndIt = remoteVerticesDeriv[rIdx+i].end();
                for (auto cIt=remoteVerticesDeriv[rIdx+i].begin(); cIt != cEndIt; cIt++)
                    cIt->umtv(nTimesEdge[i],linIntElem[0][cIt.index()][0]);
            }
            ctype intElem = (dim-1)*rGlobGeom.volume();
    
            // store integral of the products between mortar and dual functions
            std::vector<VectorType> sfIntM(nmSize),sfIntNM(nmSize);
            RangeType gapValues(nmSize);
            gapValues = 0;
            // gap vector contains the gap times mortar basis function times dual function
            std::vector<VectorType> gapVectors(nmSize);
            for (size_t k=0; k<nmSize; k++) {
                sfIntM[k].resize(mSize);
                sfIntNM[k].resize(nmSize);
                sfIntM[k] = 0;
                sfIntNM[k] = 0;
                gapVectors[k].resize(mSize);
                gapVectors[k] = 0;
            }
    
            // loop over quadrature points
            for (size_t l=0; l<quad.size(); l++) {
    
                // (dim-1)-quadrature point in the dim-reference element coordinates
                auto nonmortarQuadPos = rGeomIn.global(quad[l].position());
                auto mortarQuadPos = rGeomOut.global(quad[l].position());
    
                // QP on the nm face
                auto nmFaceQP = nmLocalFaceGeometry.local(nonmortarQuadPos);
                auto mFaceQP = mLocalFaceGeometry.local(mortarQuadPos);
    
                //evaluate shape functions at the quadrature point
                std::vector<FieldMatrix<field_type, 1, dim-1> > nmRefFaceJacobians, mRefFaceJacobians, dualRefFaceJacobians;
                std::vector<FieldVector<field_type, 1> > mFaceValues, nmFaceValues, dualFaceValues;
    
                remoteLfe.evaluateFunction(quad[l].position(),remoteQuadValues);
                mortarLfe.evaluateFunction(mortarQuadPos,mQuadValues);
                mortarLfe.evaluateJacobian(mortarQuadPos,mRefJacobians);
                nonmortarLfe.evaluateFunction(nonmortarQuadPos,nmQuadValues);
                nonmortarLfe.evaluateJacobian(nonmortarQuadPos,nmRefJacobians);
                dualLfe.evaluateFunction(nonmortarQuadPos,dualQuadValues);
                dualLfe.evaluateJacobian(nonmortarQuadPos,dualQuadRefJacobians);
                mFaceLfe.evaluateFunction(mFaceQP, mFaceValues);
                mFaceLfe.evaluateJacobian(mFaceQP, mRefFaceJacobians);
                nmFaceLfe.evaluateFunction(nmFaceQP, nmFaceValues);
                nmFaceLfe.evaluateJacobian(nmFaceQP, nmRefFaceJacobians);
                dualFaceLfe.evaluateFunction(nmFaceQP, dualFaceValues);
                dualFaceLfe.evaluateJacobian(nmFaceQP, dualRefFaceJacobians);
    
    
                // resort local jacobians to match the orientation within the undelying element
                size_t count(0);
                std::vector<FieldVector<field_type, dim-1> > nmFaceJacobians(nmRefFaceJacobians.size()), mFaceJacobians(mRefFaceJacobians.size());
                std::vector<FieldVector<field_type, dim-1> > dualFaceJacobians(dualRefFaceJacobians.size());
                for (size_t i=0; i<nmSize; i++) {
                  for (size_t j=0; j<nmSize; j++)
                    //if (std::fabs(nmValues[nmIdx[i]][0]-nmFaceValues[j])<1e-10) {
                    if ((nmCorners[i]-nmIntersection.geometry().corner(j)).two_norm()<1e-10) {
                      nmFaceJacobians[i] = nmRefFaceJacobians[j][0];
                      dualFaceJacobians[i] = dualRefFaceJacobians[j][0];
                      count++;
                      break;
                    }
                }
                if (count != nmFaceJacobians.size())
                  DUNE_THROW(Dune::Exception,"Didn't find all non-mortar face jacobians");
    
                count= 0;
                std::vector<int> map(mSize);
                for (size_t i=0; i<mSize; i++) {
                  for (size_t j=0; j<mSize; j++)
                    //if (std::fabs(mValues[mIdx[i]][0]-mFaceValues[j])<1e-8) {
                    if ((mCorners[i]-mIntersection.geometry().corner(j)).two_norm()<1e-10) {
                      mFaceJacobians[i] = mRefFaceJacobians[j][0];
                      map[i] = j;
                      count++;
                      break;
                    }
                }
    
                if (count != mFaceJacobians.size())
                  DUNE_THROW(Dune::Exception,"Didn't find all mortar face jacobians");
    
                CoordinateType normalQp;
                mortarNormal.evaluateLocal(outside,mortarQuadPos, normalQp);
                field_type normNormalQp = normalQp.two_norm();
                auto unitNormalQp = normalQp;
                unitNormalQp /= normNormalQp;
    
                // transform jacobians
                CoordinateVectorType nmJacobians(dim), mJacobians(dim), dualJacobians(dim);
                for (size_t i=0; i<nmSize;i++) {
    
                    const auto& invJacobian = nmElemGeom.jacobianInverseTransposed(nonmortarQuadPos);
    
                    invJacobian.mv(nmRefJacobians[nmIdx[i]][0], nmJacobians[i]);
                    invJacobian.mv(dualQuadRefJacobians[nmIdx[i]][0], dualJacobians[i]);
                }
    
                for (size_t i=0; i<mSize;i++) {
                    const auto& invJacobian = mElemGeom.jacobianInverseTransposed(mortarQuadPos);
                    invJacobian.mv(mRefJacobians[mIdx[i]][0], mJacobians[i]);
                }
    
                for (size_t j=0; j<nmSize; j++) {
    
                    // the integrals of the products of dual and nonmortar functions
                    for (size_t k=0; k<nmSize; k++)
                        sfIntNM[j][k].axpy(intElem*dualQuadValues[nmIdx[j]]*nmQuadValues[nmIdx[k]]*quad[l].weight(),unitNormalQp);
    
                    // the integrals of the products of dual and mortar functions
                    for (size_t k=0; k<mSize; k++)
                        sfIntM[j][k].axpy(intElem*dualQuadValues[nmIdx[j]]*mQuadValues[mIdx[k]]*quad[l].weight(),unitNormalQp);
                }
    
                // The gap function at quadrature point
                auto gapVector = nmElemGeom.global(nonmortarQuadPos) - mElemGeom.global(mortarQuadPos);
    
                for (size_t j=0; j<nmSize; j++) {
                  for (size_t k=0; k<mSize; k++)
                    gapVectors[j][k].axpy(quad[l].weight()*intElem*dualQuadValues[nmIdx[j]]*mQuadValues[mIdx[k]]/normNormalQp,gapVector);
                  gapValues[j] += quad[l].weight()*dualQuadValues[nmIdx[j]]*(gapVector*unitNormalQp);
                }
    
                ////////////////////////////////////////////////////////////
                ///          Derivation of local geometries Part A
                ///////////////////////////////////////////////////////////////
    
                /////////////////////////////////////////////////////////
                ///     To compute the linearisation of the geometryInInside method
                ///     we derive the following equation:
                ///     nonmortarGlobalQP = rGlobQP
                /// //////////////////////////////////////////////////////////
    
                // the index pattern of the non-mortar linearised geometry
                // and the integration element are the same
                MatrixType linNmGeometry;
                Dune::MatrixIndexSet nmGeomIndexSet(1,size0+size1);
                for (int i=0; i<dim; i++) {
                    const auto& row = remoteVerticesDeriv[rIdx+i];
                    const auto& cEndIt = row.end();
                    for (auto cIt=row.begin(); cIt != cEndIt; cIt++)
                        nmGeomIndexSet.add(0,cIt.index());
                }
    
                for (size_t i=0; i<nmSize; i++)
                    nmGeomIndexSet.add(0,nmglobalIdx[i]);
    
                nmGeomIndexSet.exportIdx(linNmGeometry);
                linNmGeometry = 0;
    
                // from the product rule of the left side
                for (size_t i=0; i<nmSize; i++)
                    MV::addToDiagonal(linNmGeometry[0][nmglobalIdx[i]],-nmQuadValues[nmIdx[i]]);
    
                for (int i=0; i<dim; i++) {
    
                    const auto& cEndIt = remoteVerticesDeriv[rIdx+i].end();
                    for (auto cIt=remoteVerticesDeriv[rIdx+i].begin(); cIt != cEndIt; cIt++)
                        linNmGeometry[0][cIt.index()].axpy(remoteQuadValues[i],*cIt);
                }
    
                // compute dual of covariant nonmortar basis
                auto nmDualCovBasis = dualCovBasis(nmFaceJacobians,nmCorners);
    
                //////////////////////////////////////////////////////
                /// Local Mortar: derive proj. MortarQP = remote QP
                ////////////////////////////////////////////////////
                Dune::MatrixIndexSet mGeomIndexSet(1,size0+size1);
                mGeomIndexSet.import(linNmGeometry);
                for (size_t i=0; i<mSize; i++) {
                    const auto& cEndIt = linVertex[i][0].end();
                    for (auto cIt=linVertex[i][0].begin(); cIt != cEndIt; cIt++)
                        mGeomIndexSet.add(0,cIt.index());
                }
    
                MatrixType linMGeometry;
                mGeomIndexSet.exportIdx(linMGeometry);
                linMGeometry = 0;
    
    
                // from the product rule of the left side
                for (size_t i=0; i<mSize; i++) {
    
                    const auto& cEnd = linVertex[i][0].end();
                    for (auto cI = linVertex[i][0].begin(); cI != cEnd; cI++)
                        linMGeometry[0][cI.index()].axpy(-mQuadValues[mIdx[i]],*cI);
                }
    
                for (int i=0; i<dim; i++) {
    
                    const auto& cEndIt = remoteVerticesDeriv[rIdx+i].end();
                    for (auto cIt=remoteVerticesDeriv[rIdx+i].begin(); cIt != cEndIt; cIt++)
                        linMGeometry[0][cIt.index()].axpy(remoteQuadValues[i],*cIt);
                }
    
                auto mDualCovBasis = dualCovBasis(mFaceJacobians, projVertex);
    
                /////////////////////////////////////////////////////////////////
                // Derivation of local geometries Part B
                /////////////////////////////////////////////////////////////////////////
    
                // loop over the dual functions
                for (size_t j=0; j<nmSize; j++) {
    
                    // get global index of the dual shape function
                    int rowIdx = globalToLocal[nmglobalIdx[j]];
    
                    /////////////////////////////////////////////////////////////////////////
                    // part from the local nonmortar geometry, i.e. the nonmortar matrix part
                    ///////////////////////////////////////////////////////////////////////
    
                    const auto& cNmEndIt = linNmGeometry[0].end();
                    for (auto cIt = linNmGeometry[0].begin();cIt != cNmEndIt; cIt++) {
    
                        CoordinateType vecNm(0),vecDual(0);
                        for (size_t p=0; p<nmDualCovBasis.size(); p++) {
                            CoordinateType dummy(0);
                            cIt->mtv(nmDualCovBasis[p],dummy);
                            vecDual.axpy(dualFaceJacobians[j][p],dummy);
                            vecNm.axpy(nmFaceJacobians[j][p],dummy);
                        }
    
                        for (size_t k=0; k<nmSize; k++) {
    
                            int rowIdxk = globalToLocal[nmglobalIdx[k]];
    
                            // non-mortar corner deformation times the corresponding averaged normal
                            ctype scal = (unitNormalQp*nmCorners[k])*quad[l].weight()*intElem;
    
                            linearMat[rowIdx][cIt.index()][0].axpy(scal*nmQuadValues[nmIdx[k]],vecDual);
    
                            scal = (unitNormalQp*nmCorners[j])*quad[l].weight()*intElem;
    
                            // Note: Only because the dual and lagrange functions are orthogonal,
                            // doesn't mean that also their linearisation is diagonal!
                            linearMat[rowIdxk][cIt.index()][0].axpy(scal*dualQuadValues[nmIdx[k]],vecNm);
                        }
    
                        //////////////////////////////////////////////////////
                        //          parts from the local mortar geometry
                        //////////////////////////////////////////////////////////////
    
                        //   Derivation of the dual function
                        // loop over the mortar vertices
                        for (size_t k=0; k<mSize; k++) {
    
                            // mortar corner deformation times the corresponding averaged normal
                            ctype scal = -1*(unitNormalQp*mCorners[k])*quad[l].weight()*intElem;
    
                            linearMat[rowIdx][cIt.index()][0].axpy(scal*mQuadValues[mIdx[k]],vecDual);
                        }
                    }
    
                    ////////////////////////////////////////////
                    //  Derivation of the mortar geometry
                    ////////////////////////////////////////////
    
                    const auto& cEndMGeom = linMGeometry[0].end();
                    for (size_t k=0; k<mSize; k++) {
    
                        ctype scal = -(unitNormalQp*mCorners[k])*quad[l].weight()*intElem*dualQuadValues[nmIdx[j]];
                        ctype scalNorm =(gapVector*avNormalsM[mglobalIdx[k]-size0])*quad[l].weight()*intElem*dualQuadValues[nmIdx[j]]/normNormalQp;
    
                        auto scalNor = -(normalQp*avNormalsM[mglobalIdx[k]-size0])*quad[l].weight()*intElem*dualQuadValues[nmIdx[j]]/std::pow(normNormalQp,3);
                        for (auto cI = linMGeometry[0].begin(); cI != cEndMGeom;cI++) {
    
                            CoordinateType vec(0);
                            for (size_t p=0; p<mDualCovBasis.size(); p++) {
                                CoordinateType dummy(0);
                                cI->mtv(mDualCovBasis[p],dummy);
                                vec.axpy(mFaceJacobians[k][p],dummy);
                            }
                            linearMat[rowIdx][cI.index()][0].axpy(scal,vec);
    
                            // the normal derivative
                            linearMat[rowIdx][cI.index()][0].axpy(scalNorm,vec);
    
                            // first part of the linearisation of the averaged normal field
                            // at the quadrature point the normal has to be norm again and
                            // this division has to be linearised to
                            vec *= scalNor;
                            auto linNorm = dyadic(normalQp,vec);
                            linNorm.umtv(gapVector,linearMat[rowIdx][cI.index()][0]);
                        }
    
                        // second part of linearisation of the average normal norm
                        // derivation of norm includes the average normal
                        int normalIdx = globalToLocalM[mglobalIdx[k]-size0];
                        const auto& cEnd = avNormMortarDeriv[normalIdx].end();
                        // The linearised normals are computed w.r.t local indices
                        for (auto cI = avNormMortarDeriv[normalIdx].begin(); cI != cEnd; cI++) {
                            CoordinateType avNorLin;
                            cI->mtv(normalQp, avNorLin);
                            avNorLin *= -mQuadValues[mIdx[k]]*dualQuadValues[nmIdx[j]]*quad[l].weight()*intElem/std::pow(normNormalQp,3);
    
                            auto linNorm = dyadic(normalQp,avNorLin);
                            linNorm.umtv(gapVector, linearMat[rowIdx][size0 + cI.index()][0]);
                        }
                    }
                }
            }
    
            /////////////////////////////////////////////
            //  Derivation of the av. normals part
            //////////////////////////////////////////////
    
            for (size_t j=0; j<nmSize; j++)
                for (size_t k=0; k<mSize; k++) {
    
                // get global index of the shape function
                int rowIdx = globalToLocal[nmglobalIdx[j]];
                int normalIdx = globalToLocalM[mglobalIdx[k]-size0];
    
                const auto& cEnd = avNormMortarDeriv[normalIdx].end();
                // The linearised normals are computed w.r.t local indices
                for (auto cI = avNormMortarDeriv[normalIdx].begin(); cI != cEnd; cI++)
                    cI->umtv(gapVectors[j][k],linearMat[rowIdx][size0 + cI.index()][0]);
            }
    
            ///////////////////////////////////////////////////////////////////////////////
            //      Derivation of the integration element and the deformation coefficients
            /////////////////////////////////////////////////////////////////////////////////
    
            const auto& cEnd = linIntElem[0].end();
    
            for (size_t j=0; j<nmSize; j++) {
                int rowIdx = globalToLocal[nmglobalIdx[j]];
                for (auto cIt = linIntElem[0].begin(); cIt!= cEnd; cIt++)
                  linearMat[rowIdx][cIt.index()][0].axpy(gapValues[j],(*cIt)[0]);
            }
    
            // loop over all nonmortar shape functions
            for (size_t j=0; j<nmSize; j++) {
    
                // get global index of the shape function
                int rowIdx = globalToLocal[nmglobalIdx[j]];
    
                // derivation w.r.t. the deformation  coefficient,
                //  minus because we substract the nonmortar part
                for (size_t k=0; k<nmSize; k++)
                    linearMat[rowIdx][nmglobalIdx[k]][0] += sfIntNM[j][k];
    
                // mortar part
                for (size_t k=0; k<mSize; k++)
                    linearMat[rowIdx][mglobalIdx[k]][0] -= sfIntM[j][k];
            }
        }
        std::cout<<"Computed exact derivative of mortar constraint!\n";
    }
    
    } /* namespace Contact */
    } /* namespace Dune */