Select Git revision
gapfunction.cc
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 */