diff --git a/dune/tectonic/factories/twoblocksfactory.cc b/dune/tectonic/factories/twoblocksfactory.cc
index f58defd7d1f47874e282569ec9b8ddc228173cf3..39c47cd7ab87bb3721d74ef8d84417628d59f2d6 100644
--- a/dune/tectonic/factories/twoblocksfactory.cc
+++ b/dune/tectonic/factories/twoblocksfactory.cc
@@ -35,8 +35,6 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1], depths[1]);
cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
- cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1] + heights[1], 0});
-
#elif MY_DIM == 2
origins[0] = {0, 0};
origins[1] = {lengths[0]/2.0, origins[0][1] + heights[0]};
@@ -45,8 +43,7 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
cuboidGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0]}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0]});
cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1]);
- cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1]});
- cuboidGeometries_[1]->addWeakeningPatch(frictionParset, {origins[1][0] + lengths[1], origins[1][1]}, {origins[1][0] + lengths[1], origins[1][1] + heights[1]});
+ cuboidGeometries_[1]->addWeakeningPatch(frictionParset, {origins[1][0] + 0.2*lengths[1], origins[1][1]}, {origins[1][0] + 0.8*lengths[1], origins[1][1]});
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
@@ -55,13 +52,15 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
auto& grids = gridConstructor_->getGrids();
+ std::array<double, 2> smallestDiameter = {this->parset_.template get<double>("body0.smallestDiameter"), this->parset_.template get<double>("body1.smallestDiameter")};
+
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
// define weak patch and refine grid
const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
for (size_t j=0; j<weakeningRegions.size(); j++) {
- refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::lengthScale());
+ refine(*grids[i], weakeningRegions[j], smallestDiameter[i], CuboidGeometry::lengthScale());
}
// determine minDiameter and maxDiameter
@@ -177,7 +176,7 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
// body1: Neumann boundary (upper)
// normal load
- std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+ /*std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
for (int j=0; j<leafVertexCount1; j++) {
if (leafFaces_[1]->upper.containsVertex(j))
(*loadNeumannNodes)[j][1] = true;
@@ -193,7 +192,7 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
loadNeumannBoundary->setBoundaryFunction(constantFunction);
- body1->addBoundaryCondition(loadNeumannBoundary);
+ body1->addBoundaryCondition(loadNeumannBoundary);*/
}
#include "twoblocksfactory_tmpl.cc"
diff --git a/dune/tectonic/spatial-solving/CMakeLists.txt b/dune/tectonic/spatial-solving/CMakeLists.txt
index 75fd6237669e0022c8304dc463775b465c0887b7..e7b68c8644354e5aa46db465b9dd7ef895df2e4c 100644
--- a/dune/tectonic/spatial-solving/CMakeLists.txt
+++ b/dune/tectonic/spatial-solving/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory("contact")
add_subdirectory("tnnmg")
add_subdirectory("preconditioners")
diff --git a/dune/tectonic/spatial-solving/contact/CMakeLists.txt b/dune/tectonic/spatial-solving/contact/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b651632d706f6b55dee76d0cd346e8f7e353238c
--- /dev/null
+++ b/dune/tectonic/spatial-solving/contact/CMakeLists.txt
@@ -0,0 +1,15 @@
+add_custom_target(tectonic_dune_spatial-solving_contact SOURCES
+ dualmortarcoupling.hh
+ dualmortarcoupling.cc
+ nbodyassembler.hh
+ nbodyassembler.cc
+ nbodycontacttransfer.hh
+ nbodycontacttransfer.cc
+)
+
+#install headers
+install(FILES
+ dualmortarcoupling.hh
+ nbodyassembler.hh
+ nbodycontacttransfer.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9f4daa82f20240a3232672eacee18d587555175b
--- /dev/null
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
@@ -0,0 +1,344 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#include <dune/common/exceptions.hh>
+
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+
+#include <dune/grid-glue/merging/contactmerge.hh>
+//#include <dune/grid-glue/adapter/gridgluevtkwriter.hh>
+//#include <dune/grid-glue/adapter/gridglueamirawriter.hh>
+
+#include <dune/contact/common/dualbasisadapter.hh>
+#include <dune/matrix-vector/addtodiagonal.hh>
+
+namespace Dune {
+namespace Contact {
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagrangeMatrix()
+{
+ // Create mapping from the global set of block dofs to the ones on the contact boundary
+ std::vector<int> globalToLocal;
+ nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+ // clear matrix
+ nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nonmortarBoundary_.numVertices());
+ nonmortarLagrangeMatrix_ = 0;
+
+ const auto& indexSet = grid0_->leafIndexSet();
+
+ // loop over all faces of the boundary patch
+ for (const auto& nIt : nonmortarBoundary_) {
+
+ const auto& geometry = nIt.geometry();
+ const field_type numCorners = geometry.corners();
+
+ ctype intElem = geometry.integrationElement(Dune::FieldVector<ctype,dim-1>(0));
+
+ field_type sfI = (numCorners==3) ? intElem/6.0 : intElem/numCorners;
+
+ // turn scalar element mass matrix into vector-valued one
+ // and add element mass matrix to the global matrix
+
+ // Get global node ids
+ const auto& inside = nIt.inside();
+ const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+
+ // Add localMatrix to nonmortarLagrangeMat
+ for (int i=0; i<refElement.size(nIt.indexInInside(),1,dim); i++) {
+ // we can use subEntity here because we add all indices anyway
+ int v = globalToLocal[indexSet.subIndex(inside,refElement.subEntity(nIt.indexInInside(),1,
+ i,dim),dim)];
+ Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
+ }
+
+ }
+
+}
+
+template <class field_type, class GridType0, class GridType1>
+void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
+{
+ typedef typename GridType0::LeafGridView GridView0;
+ typedef typename GridType1::LeafGridView GridView1;
+
+ typedef Dune::PQkLocalFiniteElementCache<typename GridType1::ctype, field_type, GridType1::dimension,1> FiniteElementCache1;
+
+ // cache for the dual functions on the boundary
+ using DualCache = Dune::Contact::DualBasisAdapter<GridView0, field_type>;
+
+ using Element0 = typename GridView0::template Codim<0>::Entity;
+ using Element1 = typename GridView1::template Codim<0>::Entity;
+
+ auto desc0 = [&] (const Element0& e, unsigned int face) {
+ return nonmortarBoundary_.contains(e,face);
+ };
+
+ auto desc1 = [&] (const Element1& e, unsigned int face) {
+ return mortarBoundary_.contains(e,face);
+ };
+
+ GridView0 gridView0 = grid0_->leafGridView();
+ GridView1 gridView1 = grid1_->leafGridView();
+
+ auto extract0 = std::make_shared<Extractor0>(gridView0,desc0);
+ auto extract1 = std::make_shared<Extractor1>(gridView1,desc1);
+
+ if (!gridGlueBackend_)
+ gridGlueBackend_ = std::make_shared< Dune::GridGlue::ContactMerge<dimworld, ctype> >(overlap_);
+ glue_ = std::make_shared<Glue>(extract0, extract1, gridGlueBackend_);
+ auto& glue = *glue_;
+ glue.build();
+
+ //std::cout << glue.size() << " remote intersections found." << std::endl;
+ //GridGlueAmiraWriter::write<GlueType>(glue,debugPath_);
+
+ // Restrict the hasObstacle fields to the part of the nonmortar boundary
+ // where we could actually create a contact mapping
+ LeafBoundaryPatch0 boundaryWithMapping(gridView0);
+
+ const auto& indexSet0 = gridView0.indexSet();
+ const auto& indexSet1 = gridView1.indexSet();
+
+ ///////////////////////////////////
+ // reducing nonmortar boundary
+ /////////////////////////////////
+
+ // Get all fine grid boundary segments that are totally covered by the grid-glue segments
+ typedef std::pair<int,int> Pair;
+ std::map<Pair,ctype> coveredArea, fullArea;
+
+ // initialize with area of boundary faces
+ for (const auto& bIt : nonmortarBoundary_) {
+ const Pair p(indexSet0.index(bIt.inside()),bIt.indexInInside());
+ fullArea[p] = bIt.geometry().volume();
+ coveredArea[p] = 0;
+ }
+
+ // sum up the remote intersection areas to find out which are totally covered
+ for (const auto& rIs : intersections(glue))
+ coveredArea[Pair(indexSet0.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 : nonmortarBoundary_) {
+ const auto& inside = bIt.inside();
+ if(coveredArea[Pair(indexSet0.index(inside),bIt.indexInInside())]/
+ fullArea[Pair(indexSet0.index(inside),bIt.indexInInside())] >= coveredArea_)
+ boundaryWithMapping.addFace(inside, bIt.indexInInside());
+ }
+
+ //writeBoundary(boundaryWithMapping,debugPath_ + "relevantNonmortar");
+
+
+ /** \todo replace by all fine grid segments which are totally covered by the RemoteIntersections. */
+ //for (const auto& rIs : intersections(glue))
+ // boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());
+
+ /*printf("contact mapping could be built for %d of %d boundary segments.\n",
+ boundaryWithMapping.numFaces(), nonmortarBoundary_.numFaces()); */
+
+ nonmortarBoundary_ = boundaryWithMapping;
+ mortarBoundary_.setup(gridView1);
+ for (const auto& rIs : intersections(glue))
+ if (nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ mortarBoundary_.addFace(rIs.outside(),rIs.indexInOutside());
+
+
+ // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
+ assembleNonmortarLagrangeMatrix();
+
+ // The weak obstacle vector
+ weakObstacle_.resize(nonmortarBoundary_.numVertices());
+ weakObstacle_ = 0;
+
+ // ///////////////////////////////////////////////////////////
+ // Get the occupation structure for the mortar matrix
+ // ///////////////////////////////////////////////////////////
+
+ /** \todo Also restrict mortar indices and don't use the whole grid level. */
+ Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));
+
+ // Create mapping from the global set of block dofs to the ones on the contact boundary
+ std::vector<int> globalToLocal;
+ nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+ // loop over all intersections
+ for (const auto& rIs : intersections(glue)) {
+
+ if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ continue;
+
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+ const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+ int nDomainVertices = domainRefElement.size(dim);
+ int nTargetVertices = targetRefElement.size(dim);
+
+ for (int j=0; j<nDomainVertices; j++) {
+
+ int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];
+
+ // if the vertex is not contained in the restricted contact boundary then dismiss it
+ if (localDomainIdx == -1)
+ continue;
+
+ for (int k=0; k<nTargetVertices; k++) {
+ int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
+ if (!mortarBoundary_.containsVertex(globalTargetIdx))
+ continue;
+
+ mortarIndices.add(localDomainIdx, globalTargetIdx);
+ }
+ }
+ }
+
+ mortarIndices.exportIdx(mortarLagrangeMatrix_);
+
+ // Clear it
+ mortarLagrangeMatrix_ = 0;
+
+
+ //cache of local bases
+ FiniteElementCache1 cache1;
+
+ std::unique_ptr<DualCache> dualCache;
+ dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView0, field_type> >();
+
+ std::vector<Dune::FieldVector<ctype,dim> > avNormals;
+ avNormals = nonmortarBoundary_.getNormals();
+
+ // loop over all intersections and compute the matrix entries
+ for (const auto& rIs : intersections(glue)) {
+
+
+ const auto& inside = rIs.inside();
+
+ if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ continue;
+
+ const auto& outside = rIs.outside();
+
+ // types of the elements supporting the boundary segments in question
+ Dune::GeometryType nonmortarEType = inside.type();
+ Dune::GeometryType mortarEType = outside.type();
+
+ const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+
+ const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(mortarEType);
+
+ int noOfMortarVec = targetRefElement.size(dim);
+
+ 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).
+ int quadOrder = 2 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
+ const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(rIs.type(), quadOrder);
+
+ //const typename FiniteElementCache0::FiniteElementType& nonmortarFiniteElement = cache0.get(nonmortarEType);
+ const auto& mortarFiniteElement = cache1.get(mortarEType);
+ dualCache->bind(inside, rIs.indexInInside());
+
+ std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;
+
+ const auto& rGeom = rIs.geometry();
+ const auto& rGeomOutside = rIs.geometryOutside();
+ const auto& rGeomInInside = rIs.geometryInInside();
+ const auto& rGeomInOutside = rIs.geometryInOutside();
+
+ int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
+ std::vector<int> nonmortarFaceNodes;
+ for (int i=0; i<nNonmortarFaceNodes; i++) {
+ int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
+ nonmortarFaceNodes.push_back(faceIdxi);
+ }
+
+ for (const auto& quadPt : quadRule) {
+
+ // compute integration element of overlap
+ ctype integrationElement = rGeom.integrationElement(quadPt.position());
+
+ // quadrature point positions on the reference element
+ Dune::FieldVector<ctype,dim> nonmortarQuadPos = rGeomInInside.global(quadPt.position());
+ Dune::FieldVector<ctype,dim> mortarQuadPos = rGeomInOutside.global(quadPt.position());
+
+ // The current quadrature point in world coordinates
+ Dune::FieldVector<field_type,dim> nonmortarQpWorld = rGeom.global(quadPt.position());
+ Dune::FieldVector<field_type,dim> mortarQpWorld = rGeomOutside.global(quadPt.position());;
+
+ // the gap direction (normal * gapValue)
+ Dune::FieldVector<field_type,dim> gapVector = mortarQpWorld - nonmortarQpWorld;
+
+ //evaluate all shapefunctions at the quadrature point
+ //nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
+ mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
+ dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+ // loop over all Lagrange multiplier shape functions
+ for (int j=0; j<nNonmortarFaceNodes; j++) {
+
+ int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+ int rowIdx = globalToLocal[globalDomainIdx];
+
+ weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+ * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+
+ // loop over all mortar shape functions
+ for (int k=0; k<noOfMortarVec; k++) {
+
+ int colIdx = indexSet1.subIndex(outside, k, dim);
+ if (!mortarBoundary_.containsVertex(colIdx))
+ continue;
+
+ // Integrate over the product of two shape functions
+ field_type mortarEntry = integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+
+ Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+ }
+
+ }
+
+ }
+
+ }
+
+ // ///////////////////////////////////////
+ // Compute M = D^{-1} \hat{M}
+ // ///////////////////////////////////////
+
+ Dune::BCRSMatrix<MatrixBlock>& M = mortarLagrangeMatrix_;
+ Dune::BDMatrix<MatrixBlock>& D = nonmortarLagrangeMatrix_;
+
+ // First compute D^{-1}
+ D.invert();
+
+ // Then the matrix product D^{-1} \hat{M}
+ for (auto rowIt = M.begin(); rowIt != M.end(); ++rowIt) {
+ const auto rowIndex = rowIt.index();
+ for (auto& entry : *rowIt)
+ entry.leftmultiply(D[rowIndex][rowIndex]);
+ }
+
+ // weakObstacles in transformed basis = D^{-1}*weakObstacle_
+ for(size_t rowIdx=0; rowIdx<weakObstacle_.size(); rowIdx++)
+ weakObstacle_[rowIdx] *= D[rowIdx][rowIdx][0][0];
+
+ gridGlueBackend_->clear();
+}
+
+} /* namespace Contact */
+} /* namespace Dune */
diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cfd8ed5f25e46034afadc2acdbccd56f386a8e48
--- /dev/null
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh
@@ -0,0 +1,192 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_CONTACT_ASSEMBLERS_DUAL_MORTAR_COUPLING_HH
+#define DUNE_CONTACT_ASSEMBLERS_DUAL_MORTAR_COUPLING_HH
+
+#include <memory>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bdmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/boundarypatchprolongator.hh>
+
+#include <dune/grid-glue/gridglue.hh>
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+#include <dune/grid-glue/merging/merger.hh>
+
+namespace Dune {
+namespace Contact {
+
+/** \brief Assembles the transfer operator for two-body contact
+ */
+template<class field_type, class GridType0, class GridType1=GridType0>
+class DualMortarCoupling {
+
+ // /////////////////////
+ // private types
+ // /////////////////////
+
+ enum {dim = GridType0::dimension};
+
+ enum {dimworld = GridType0::dimensionworld};
+
+ protected:
+ typedef Dune::FieldMatrix<field_type, dim, dim> MatrixBlock;
+ typedef Dune::BCRSMatrix<MatrixBlock> MatrixType;
+ using ObstacleVector = Dune::BlockVector<Dune::FieldVector<field_type, 1> >;
+
+ typedef typename GridType0::ctype ctype;
+
+ typedef BoundaryPatch<typename GridType0::LeafGridView> LeafBoundaryPatch0;
+ typedef BoundaryPatch<typename GridType1::LeafGridView> LeafBoundaryPatch1;
+
+ typedef BoundaryPatch<typename GridType0::LevelGridView> LevelBoundaryPatch0;
+ typedef BoundaryPatch<typename GridType1::LevelGridView> LevelBoundaryPatch1;
+
+ using Extractor0 = Dune::GridGlue::Codim1Extractor<typename GridType0::LeafGridView>;
+ using Extractor1 = Dune::GridGlue::Codim1Extractor<typename GridType1::LeafGridView>;
+
+public:
+
+ using Glue = Dune::GridGlue::GridGlue<Extractor0, Extractor1>;
+
+ DualMortarCoupling(field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
+ : overlap_(overlap), gridGlueBackend_(nullptr), coveredArea_(coveredArea)
+ {}
+
+ DualMortarCoupling(const GridType0& grid0, const GridType1& grid1,
+ field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
+ : grid0_(&grid0), grid1_(&grid1), overlap_(overlap),
+ gridGlueBackend_(nullptr), coveredArea_(coveredArea)
+ {}
+
+ virtual ~DualMortarCoupling() {
+ /* Nothing. */
+ }
+
+ /** \brief Sets up the contact coupling operator on the grid leaf level */
+ virtual void setup();
+
+ /** \brief Assemble nonmortar matrix on the leaf level. */
+ void assembleNonmortarLagrangeMatrix();
+
+ /** \brief Get the dual mortar coupling matrix of the leaf level. */
+ virtual const MatrixType& mortarLagrangeMatrix() const {
+ return mortarLagrangeMatrix_;
+ }
+
+ /** \brief Get the dual mortar coupling matrix of the leaf level. */
+ virtual MatrixType& mortarLagrangeMatrix() {
+ return mortarLagrangeMatrix_;
+ }
+
+ /** \brief Get the non-mortar matrix. */
+ const BDMatrix<MatrixBlock>& nonmortarLagrangeMatrix() const {
+ return nonmortarLagrangeMatrix_;
+ }
+
+ /** \brief Setup leaf nonmortar and mortar patches. */
+ virtual void setupContactPatch(const LevelBoundaryPatch0& coarseNonmortar,
+ const LevelBoundaryPatch1& coarseMortar) {
+
+ BoundaryPatchProlongator<GridType0>::prolong(coarseNonmortar,nonmortarBoundary_);
+ BoundaryPatchProlongator<GridType1>::prolong(coarseMortar,mortarBoundary_);
+ }
+
+ /** \brief Return the leafnonmortar boundary. */
+ const LeafBoundaryPatch0& nonmortarBoundary() const {
+ return nonmortarBoundary_;
+ }
+
+ /** \brief Return the leafmortar boundary. */
+ const LeafBoundaryPatch1& mortarBoundary() const {
+ return mortarBoundary_;
+ }
+
+ /** \brief Set the non-mortar and mortar grid. */
+ void setGrids(const GridType0& grid0, const GridType1& grid1) {
+ grid0_ = &grid0;
+ grid1_ = &grid1;
+ }
+
+ /** \brief Get non-mortar grid. */
+ const GridType0* getGrid0() const {
+ return grid0_;
+ }
+
+ /** \brief Get mortar grid. */
+ const GridType1* getGrid1() const {
+ return grid1_;
+ }
+
+ /** \brief Set the percentage of covered area each face must at least have
+ * by the grid-glue projection.
+ */
+ void setCoveredArea(field_type coveredArea) {
+ coveredArea_ = coveredArea;
+ }
+
+ void setOverlap(field_type overlap) {
+ overlap_ = overlap;
+ }
+
+ /** \brief Get the GridGlue. */
+ std::shared_ptr<Glue> getGlue() {
+ return glue_;
+ }
+
+ const ObstacleVector& getWeakObstacle() const {
+ return weakObstacle_;
+ }
+
+ // /////////////////////////////////////////////
+ // Data members
+ // /////////////////////////////////////////////
+
+ /** \brief The two grids involved in the two-body contact problem. */
+ const GridType0* grid0_;
+ const GridType1* grid1_;
+
+ /** \brief For each dof a bit specifying whether the dof carries an obstacle or not. */
+ LeafBoundaryPatch0 nonmortarBoundary_;
+
+ /** \brief The mortar boundary. */
+ LeafBoundaryPatch1 mortarBoundary_;
+
+ /** \brief The crosspoints of the nonmortar body. */
+ std::set<size_t> crosspoints_;
+
+ /** \brief The matrix coupling mortar side and Lagrange multipliers. */
+ MatrixType mortarLagrangeMatrix_;
+
+ /** \brief The diagonal matrix coupling nonmortar side and Lagrange multipliers */
+ Dune::BDMatrix<MatrixBlock> nonmortarLagrangeMatrix_;
+
+ /** \brief The weak obstacles. */
+ ObstacleVector weakObstacle_;
+
+ /** \brief Allow some small penetration of the bodies. */
+ field_type overlap_;
+
+ /** \brief A grid-glue merger. If this is not set ContactMerge is used. */
+ std::shared_ptr< Dune::GridGlue::Merger<field_type,dim-1,dim-1,dim> > gridGlueBackend_;
+
+ // Path where to write the merged grids and the reduced contact boundaries
+ //std::string debugPath_;
+
+ std::shared_ptr<Glue> glue_;
+
+ /** \brief Dismiss all faces that are not at least covered by the grid-glue projection for this
+ * much percentage ranging between one - for total coverage and zero for taking all faces.
+ */
+ field_type coveredArea_;
+};
+
+} /* namespace Contact */
+} /* namespace Dune */
+
+#include "dualmortarcoupling.cc"
+
+#endif
diff --git a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2d75582fa4c6cd16ae67ce340f055ccfcb10dafe
--- /dev/null
+++ b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
@@ -0,0 +1,519 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#include <memory>
+
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/matrix-vector/transpose.hh>
+#include <dune/matrix-vector/blockmatrixview.hh>
+
+#include <dune/fufem/boundarypatchprolongator.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+namespace Dune {
+namespace Contact {
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
+{
+
+ // /////////////////////////////////////////////////////////////////
+ // Check if contact data is present
+ // /////////////////////////////////////////////////////////////////
+
+ for (int i=0; i<nCouplings(); i++) {
+
+
+ if (!coupling_[i].obsPatch_)
+ DUNE_THROW(Dune::Exception, "You have to supply a nonmortar patch for the " << i << "-th coupling!");
+
+ if (!coupling_[i].mortarPatch_)
+ DUNE_THROW(Dune::Exception, "You have to supply a mortar patch for the " << i << "-th coupling!");
+ }
+
+ // ////////////////////////////////////////////////////
+ // Set up Mortar element transfer operators
+ // ///////////////////////////////////////////////////
+
+ //std::cout<<"Setup mortar transfer operators\n";
+
+ for (size_t i=0; i<contactCoupling_.size(); i++) {
+ contactCoupling_[i]->setGrids(*grids_[coupling_[i].gridIdx_[0]],*grids_[coupling_[i].gridIdx_[1]]);
+ contactCoupling_[i]->setupContactPatch(*coupling_[i].patch0(),*coupling_[i].patch1());
+ contactCoupling_[i]->gridGlueBackend_ = coupling_[i].backend();
+ contactCoupling_[i]->setCoveredArea(coveredArea_);
+ contactCoupling_[i]->setup();
+ }
+
+ // setup Householder reflections
+ assembleHouseholderReflections();
+
+ // compute the mortar transformation matrix
+ computeTransformationMatrix();
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
+{
+ // //////////////////////////////////////////////////////////////////
+ // Compute local coordinate systems for the dofs with obstacles
+ // //////////////////////////////////////////////////////////////////
+
+ std::vector<std::vector<MatrixBlock> > coordinateSystems(nCouplings());
+
+ for (int i=0; i<nCouplings(); i++) {
+
+ double dist = coupling_[i].obsDistance_;
+ auto projection = coupling_[i].projection();
+
+ if (!projection)
+ DUNE_THROW(Dune::Exception, "You have to supply a contact projection for the " << i << "-th coupling!");
+
+ std::vector<Dune::FieldVector<field_type,dim> > directions;
+
+ const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+ const auto& mortarBoundary = contactCoupling_[i]->mortarBoundary();
+
+ projection->project(nonmortarBoundary, mortarBoundary,dist);
+ projection->getObstacleDirections(directions);
+
+ this->computeLocalCoordinateSystems(nonmortarBoundary,coordinateSystems[i], directions);
+ }
+
+ // ///////////////////////////////////////////////////////////////
+ // Combine the coordinate systems for each grid to one long array
+ // ///////////////////////////////////////////////////////////////
+ this->localCoordSystems_.resize(numDofs());
+
+ // initialise with identity
+ Dune::ScaledIdentityMatrix<field_type,dim> id(1);
+ for (size_t i=0; i<this->localCoordSystems_.size(); i++)
+ this->localCoordSystems_[i] = id;
+
+ for (int j=0; j<nCouplings(); j++) {
+
+ int grid0Idx = coupling_[j].gridIdx_[0];
+
+ // compute offset
+ int idx = 0;
+ for (int k=0;k<grid0Idx;k++)
+ idx += grids_[k]->size(dim);
+
+ const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
+
+ // if a body is non-mortar more then once, one has to be careful with not overwriting entries
+ for (std::size_t k=0; k<coordinateSystems[j].size(); k++)
+ if(nonmortarBoundary.containsVertex(k))
+ this->localCoordSystems_[idx+k] = coordinateSystems[j][k];
+ }
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::assembleObstacle()
+{
+
+ std::vector<std::vector<int> > globalToLocals(nCouplings());
+ for (std::size_t i = 0; i < globalToLocals.size(); ++i)
+ contactCoupling_[i]->nonmortarBoundary().makeGlobalToLocal(globalToLocals[i]);
+
+ ///////////////////////////////////////
+ // Set the obstacle values
+ ///////////////////////////////////////
+
+ totalHasObstacle_.resize(numDofs(),false);
+
+ for (int j=0; j<nCouplings(); j++) {
+
+ int grid0Idx = coupling_[j].gridIdx_[0];
+ const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
+
+ // compute offset
+ int idx = 0;
+ for (int k=0;k<grid0Idx;k++)
+ idx += grids_[k]->size(dim);
+
+ for (int k=0; k<grids_[grid0Idx]->size(dim); k++)
+ if (nonmortarBoundary.containsVertex(k))
+ totalHasObstacle_[idx+k] = true;
+ }
+
+ // //////////////////////////////////////////////////////////////////
+ // Set the obstacle distances
+ // //////////////////////////////////////////////////////////////////
+
+ // Combine the obstacle values for each grid to one long array
+ totalObstacles_.resize(numDofs());
+
+ for (size_t j=0; j<totalObstacles_.size(); j++)
+ totalObstacles_[j].clear();
+
+ // Init the obstacle values
+ for (int i=0; i<nCouplings(); i++) {
+
+ int grid0Idx = coupling_[i].gridIdx_[0];
+
+ // compute offset
+ int idx = 0;
+ for (int k=0;k<grid0Idx;k++)
+ idx += grids_[k]->size(dim);
+
+ // The grids involved in this coupling
+ const auto& indexSet = grids_[grid0Idx]->leafIndexSet();
+
+ // the obstacles are stored using local indices
+ const std::vector<int>& globalToLocal = globalToLocals[i];
+
+ // the weak obstacles
+ const auto& obs = contactCoupling_[i]->weakObstacle_;
+
+ // get strong obstacles, the projection method of the correct level has already been called
+ //std::vector<field_type> obs;
+ //coupling_[i].projection()->getObstacles(obs);
+
+ const auto& leafView = grids_[grid0Idx]->leafGridView();
+
+ for (const auto& v : vertices(leafView)) {
+
+ int vIdx = indexSet.index(v);
+ if (globalToLocal[vIdx]<0)
+ continue;
+
+ // Set the obstacle
+ switch (coupling_[i].type_) {
+
+ case CouplingPairBase::CONTACT:
+ totalObstacles_[idx+vIdx].upper(0) = obs[globalToLocal[vIdx]];
+ break;
+
+ case CouplingPairBase::STICK_SLIP:
+ totalObstacles_[idx+vIdx].lower(0) = 0;
+ totalObstacles_[idx+vIdx].upper(0) = 0;
+ break;
+
+ case CouplingPairBase::GLUE:
+ for (int j=0; j<dim; j++) {
+ totalObstacles_[idx+vIdx].lower(j) = 0;
+ totalObstacles_[idx+vIdx].upper(j) = 0;
+ }
+ break;
+
+ case CouplingPairBase::NONE:
+ break;
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Coupling type not handled yet!");
+ }
+ }
+ }
+
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::
+assembleJacobian(const std::vector<const MatrixType*>& submat,
+ MatrixType& totalMatrix) const
+{
+ // Create a block view of the grand matrix
+ Dune::MatrixVector::BlockMatrixView<MatrixType> blockView(submat);
+
+ int nRows = blockView.nRows();
+ int nCols = blockView.nCols();
+
+ // Create the index set of B \hat{A} B^T
+ Dune::MatrixIndexSet indicesBABT(nRows, nCols);
+
+ for (int i=0; i<nGrids(); i++) {
+
+ for (size_t iRow = 0; iRow<submat[i]->N(); iRow++) {
+
+ const RowType& row = (*submat[i])[iRow];
+
+ // Loop over all columns of the stiffness matrix
+ ConstColumnIterator j = row.begin();
+ ConstColumnIterator jEnd = row.end();
+
+ for (; j!=jEnd; ++j) {
+
+ ConstColumnIterator k = BT_[blockView.row(i,iRow)].begin();
+ ConstColumnIterator kEnd = BT_[blockView.row(i,iRow)].end();
+ for (; k!=kEnd; ++k) {
+
+ ConstColumnIterator l = BT_[blockView.col(i,j.index())].begin();
+ ConstColumnIterator lEnd = BT_[blockView.col(i,j.index())].end();
+
+ for (; l!=lEnd; ++l)
+ indicesBABT.add(k.index(), l.index());
+
+ }
+
+ }
+
+ }
+
+ }
+
+ // ////////////////////////////////////////////////////////////////////
+ // Multiply transformation matrix to the global stiffness matrix
+ // ////////////////////////////////////////////////////////////////////
+
+ indicesBABT.exportIdx(totalMatrix);
+ totalMatrix = 0;
+
+ for (int i=0; i<nGrids(); i++) {
+
+ for (size_t iRow = 0; iRow<submat[i]->N(); iRow++) {
+
+ const RowType& row = (*submat[i])[iRow];
+
+ // Loop over all columns of the stiffness matrix
+ ConstColumnIterator j = row.begin();
+ ConstColumnIterator jEnd = row.end();
+
+ for (; j!=jEnd; ++j) {
+
+ ConstColumnIterator k = BT_[blockView.row(i,iRow)].begin();
+ ConstColumnIterator kEnd = BT_[blockView.row(i,iRow)].end();
+
+ for (; k!=kEnd; ++k) {
+
+ ConstColumnIterator l = BT_[blockView.col(i,j.index())].begin();
+ ConstColumnIterator lEnd = BT_[blockView.col(i,j.index())].end();
+
+ for (; l!=lEnd; ++l) {
+
+ //BABT[k][l] += BT[i][k] * mat_[i][j] * BT[j][l];
+ MatrixBlock m_ij = *j;
+
+ MatrixBlock BT_ik_trans;
+ Dune::MatrixVector::transpose(*k, BT_ik_trans);
+
+ m_ij.leftmultiply(BT_ik_trans);
+ m_ij.rightmultiply(*l);
+
+ totalMatrix[k.index()][l.index()] += m_ij;
+
+ }
+
+ }
+
+ }
+
+ }
+
+ }
+}
+
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+assembleRightHandSide(const VectorTypeContainer& rhs,
+ VectorType& totalRhs) const
+{
+ // Concatenate the two rhs vectors to a large one
+ VectorType untransformedTotalRhs(BT_.M());
+
+ int idx = 0;
+
+ for (size_t i=0; i<rhs.size(); i++) {
+
+ for (size_t j=0; j<rhs[i].size(); j++)
+ untransformedTotalRhs[idx++] = rhs[i][j];
+
+ }
+
+ if ((int) BT_.M() != idx)
+ DUNE_THROW(Dune::Exception, "assembleRightHandSide: vector size and matrix size don't match!");
+
+ // Transform the basis of the ansatz space
+ totalRhs.resize(untransformedTotalRhs.size());
+ totalRhs = 0;
+ BT_.umtv(untransformedTotalRhs, totalRhs);
+}
+
+
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+postprocess(const VectorType& totalX, VectorTypeContainer& x) const
+{
+ // ///////////////////////////////////////////////////////
+ // Transform the solution vector to the nodal basis
+ // ///////////////////////////////////////////////////////
+
+ VectorType nodalBasisTotalX(totalX.size());
+ BT_.mv(totalX, nodalBasisTotalX);
+
+
+ // ///////////////////////////////////////////////////////
+ // Split the total solution vector into the parts
+ // corresponding to the grids.
+ // ///////////////////////////////////////////////////////
+
+ int idx = 0;
+ for (int i=0; i<nGrids(); i++) {
+
+ x[i].resize(grids_[i]->size(dim));
+
+ for (size_t j=0; j<x[i].size(); j++, idx++)
+ x[i][j] = nodalBasisTotalX[idx];
+ }
+}
+
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+concatenateVectors(const VectorTypeContainer& parts, VectorType& whole)
+{
+ int totalSize = 0;
+ for (size_t i=0; i<parts.size(); i++)
+ totalSize += parts[i].size();
+
+ whole.resize(totalSize);
+
+ int idx = 0;
+ for (size_t i=0; i<parts.size(); i++)
+ for (size_t j=0; j<parts[i].size(); j++)
+ whole[idx++] = parts[i][j];
+
+}
+
+// ////////////////////////////////////////////////////////////////////////////
+// Turn the initial solutions from the nodal basis to the transformed basis
+// ////////////////////////////////////////////////////////////////////////////
+template <class GridType, class VectorType>
+template <class VectorTypeContainer>
+void NBodyAssembler<GridType, VectorType>::
+nodalToTransformed(const VectorTypeContainer& x,
+ VectorType& transformedX) const
+{
+ VectorType canonicalTotalX;
+
+ concatenateVectors(x, canonicalTotalX);
+
+ // Make small cg solver
+ Dune::MatrixAdapter<MatrixType,VectorType,VectorType> op(getTransformationMatrix());
+
+ // A preconditioner
+ Dune::SeqILU<MatrixType,VectorType,VectorType> ilu0(getTransformationMatrix(),1.0);
+
+ // A cg solver for nonsymmetric matrices
+ Dune::BiCGSTABSolver<VectorType> bicgstab(op,ilu0,1E-4,100,0);
+
+ // Object storing some statistics about the solving process
+ Dune::InverseOperatorResult statistics;
+
+ // Solve!
+ transformedX = canonicalTotalX; // seems to be a good initial value
+ bicgstab.apply(transformedX, canonicalTotalX, statistics);
+}
+
+template <class GridType, class VectorType>
+void NBodyAssembler<GridType, VectorType>::
+computeTransformationMatrix()
+{
+ //std::cout<<"Setup transformation matrix...";
+ // compute offsets for the different grids
+ std::vector<int> offsets(grids_.size());
+ offsets[0] = 0;
+
+ // P1 elements are hard-wired here
+ size_t k;
+ for (k=0; k<grids_.size()-1; k++)
+ offsets[k+1] = offsets[k] + grids_[k]->size(dim);
+
+ int nRows = offsets[k] + grids_[k]->size(dim);
+ int nCols = nRows;
+
+ // /////////////////////////////////////////////////////////////////
+ // First create the index structure
+ // /////////////////////////////////////////////////////////////////
+
+ Dune::MatrixIndexSet indicesBT(nRows, nCols);
+
+ // BT_ is the identity plus some off-diagonal elements
+ for (size_t i=0; i<indicesBT.rows(); i++)
+ indicesBT.add(i,i);
+
+ std::vector<std::vector<int> > nonmortarToGlobal(nCouplings());
+ // Enter all the off-diagonal entries
+ for (int i=0; i<nCouplings(); i++) {
+
+ const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+ // If the contact mapping could not be built at all then skip this coupling
+ if (nonmortarBoundary.numVertices() == 0)
+ continue;
+
+ // The grids involved in this coupling
+ int grid0Idx = coupling_[i].gridIdx_[0];
+ int grid1Idx = coupling_[i].gridIdx_[1];
+
+ // The mapping from nonmortar vertex indices to grid indices
+ nonmortarToGlobal[i].resize(nonmortarBoundary.numVertices());
+ int idx = 0;
+ for (int j=0; j<grids_[grid0Idx]->size(dim); j++)
+ if (nonmortarBoundary.containsVertex(j))
+ nonmortarToGlobal[i][idx++] = j;
+
+ // the off-diagonal part
+ const MatrixType& M = contactCoupling_[i]->mortarLagrangeMatrix();
+ for (size_t j=0; j<M.N(); j++) {
+
+ const RowType& row = M[j];
+
+ ConstColumnIterator cIt = row.begin();
+ ConstColumnIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt)
+ indicesBT.add(offsets[grid0Idx] + nonmortarToGlobal[i][j],
+ offsets[grid1Idx] + cIt.index());
+
+ }
+
+ }
+
+ // ////////////////////////////////////////////////////////////
+ // Enter the values of the different couplings
+ // ////////////////////////////////////////////////////////////
+
+ indicesBT.exportIdx(BT_);
+ BT_ = 0;
+
+ // Enter identity part
+ for (int i=0; i<nRows; i++)
+ for (int j=0; j<dim; j++)
+ BT_[i][i][j][j] = 1;
+
+ for (int i=0; i<nCouplings(); i++) {
+
+ const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
+ if (nonmortarBoundary.numVertices() == 0)
+ continue;
+
+ // The grids involved in this coupling
+ int grid0Idx = coupling_[i].gridIdx_[0];
+ int grid1Idx = coupling_[i].gridIdx_[1];
+
+ // the off-diagonal part
+ const MatrixType& M = contactCoupling_[i]->mortarLagrangeMatrix();
+ for (size_t j=0; j<M.N(); j++) {
+
+ const RowType& row = M[j];
+
+ ConstColumnIterator cIt = row.begin();
+ ConstColumnIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt)
+ BT_[offsets[grid0Idx] + nonmortarToGlobal[i][j]][offsets[grid1Idx] +cIt.index()] = *cIt;
+
+ }
+
+ int offset = offsets[grid0Idx];
+
+ // modify non-mortar dofs and rotate them in normal and tangential coordinates
+ for (int k=0; k<grids_[grid0Idx]->size(dim); k++)
+ if(nonmortarBoundary.containsVertex(k))
+ BT_[offset + k][offset+k] =this->localCoordSystems_[offset + k];
+ }
+}
+
+} /* namespace Contact */
+} /* namespace Dune */
diff --git a/dune/tectonic/spatial-solving/contact/nbodyassembler.hh b/dune/tectonic/spatial-solving/contact/nbodyassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..55a5d6335df2de661260620094798ff8b671a973
--- /dev/null
+++ b/dune/tectonic/spatial-solving/contact/nbodyassembler.hh
@@ -0,0 +1,206 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_CONTACT_ASSEMBLERS_N_BODY_MORTAR_ASSEMBLER_HH
+#define DUNE_CONTACT_ASSEMBLERS_N_BODY_MORTAR_ASSEMBLER_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/promotiontraits.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/contact/assemblers/dualmortarcouplinghierarchy.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include <dune/contact/assemblers/contactassembler.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/contactprojection.hh>
+
+#ifdef HAVE_DUNE_GRID_GLUE
+#include <dune/grid-glue/merging/merger.hh>
+#endif
+
+namespace Dune {
+namespace Contact {
+
+/** \brief Assembler for mortar discretized contact problems with arbitrary number of bodies.
+ *
+ * \tparam GridType Type of the corresponding grids.
+ * \tparam VectorType The vector type for the displacements.
+ */
+template <class GridType, class VectorType>
+class NBodyAssembler : public ContactAssembler<GridType::dimension, typename VectorType::field_type>
+{
+protected:
+ enum {dim = GridType::dimension};
+
+ using field_type = typename PromotionTraits<typename VectorType::field_type,
+ typename GridType::ctype>::PromotedType;
+
+ using CouplingType = DualMortarCoupling<field_type, GridType>;
+ using HierarchyCouplingType = DualMortarCouplingHierarchy<field_type, GridType>;
+
+ using MatrixBlock = Dune::FieldMatrix<field_type, dim, dim>;
+ using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
+
+ using RowType = typename MatrixType::row_type;
+ using ConstColumnIterator = typename RowType::ConstIterator;
+
+ using LeafBoundaryPatch = BoundaryPatch<typename GridType::LeafGridView>;
+
+public:
+ /** \brief Construct assembler from number of couplings and grids.
+ *
+ * \param nGrids The number of involved bodies.
+ * \param nCouplings The number of couplings.
+ * \param hierarchyCoupling This boolean determines if the whole hierarchy of mortar matrices should be setup
+ */
+ NBodyAssembler(int nGrids, int nCouplings, bool hierarchyCoupling=false, field_type coveredArea = 0.99) :
+ coveredArea_(coveredArea)
+ {
+ grids_.resize(nGrids);
+
+ coupling_.resize(nCouplings);
+ contactCoupling_.resize(nCouplings);
+ // initialise the couplings with DualMortarCoupling
+ for (int i=0; i<nCouplings; i++)
+ if (hierarchyCoupling)
+ contactCoupling_[i] = std::make_shared<HierarchyCouplingType>();
+ else
+ contactCoupling_[i] = std::make_shared<CouplingType>();
+ }
+
+ /** \brief Get the number of couplings.*/
+ int nCouplings() const {return coupling_.size();}
+
+ /** \brief Get the number of involved bodies.*/
+ int nGrids() const {return grids_.size();}
+
+ /** \brief Get the total number of degrees of freedom of the leaf grids. */
+ int numDofs() const {
+ int n=0;
+ for (int i=0; i<nGrids(); i++)
+ n += grids_[i]->size(dim);
+ return n;
+ }
+
+ /** \brief Setup the patches for the contact boundary and compute the obstacles. */
+ void assembleObstacle();
+
+ /** \brief Assemble the mortar matrices and compute the basis transformation.*/
+ void assembleTransferOperator();
+
+ /** \brief Turn initial solutions from nodal basis to the transformed basis.
+ * i.e. transformedX = O*B^{-T}nodalX
+ * \param x Initial solution vector containing nodal coefficients.
+ * \param transformedX Coefficient vector for all bodies in mortar transformed coordinates.
+ */
+ template <class VectorTypeContainer>
+ void nodalToTransformed(const VectorTypeContainer& x,
+ VectorType& transformedX) const;
+
+ /** \brief Compute stiffness matrices in mortar transformed coordinates.
+ * i.e. transformedA = O*B*nodalA*B^T*O^T
+ * \param submat Stiffness matrices w.r.t the nodal finite element basis.
+ * \param totalMatrix Reference to mortar transformed stiffness matrix for all bodies.
+ */
+ void assembleJacobian(const std::vector<const MatrixType*>& submat,
+ MatrixType& totalMatrix) const;
+
+ /** \brief Compute the right hand side in mortar transformed coordinates.
+ * i.e. totalRhs = O*B*nodalRhs
+ * \param rhs Right hand side coefficients w.r.t. the nodal finite element basis.
+ * \param totalRhs Reference to coefficient vector w.r.t. the transformed basis.
+ */
+ template <class VectorTypeContainer>
+ void assembleRightHandSide(const VectorTypeContainer& rhs,
+ VectorType& totalRhs) const;
+
+ /** \brief Transform a vector from local coordinates to canonical coordinates.
+ *
+ * \param totalX Coefficient vector of the mortar transformed basis.
+ * \param x Reference to target vector for the standard nodal coefficients of each body.
+ */
+ template <class VectorTypeContainer>
+ void postprocess(const VectorType& totalX, VectorTypeContainer& x) const;
+
+ /** \brief Concatenate a family of vectors .
+ *
+ * \param parts A vector of vectors.
+ * \param whole The vector to contain the concatenated family
+ */
+ template <class VectorTypeContainer>
+ static void concatenateVectors(const VectorTypeContainer& parts, VectorType& whole);
+
+ /** \brief Get the contact couplings. */
+ const auto& getContactCouplings() const {
+ return contactCoupling_;
+ }
+
+ /** \brief Get the contact couplings. */
+ auto& getContactCouplings() {
+ return contactCoupling_;
+ }
+
+ /** \brief Set the contact couplings. */
+ void setContactCouplings(std::vector<std::shared_ptr<CouplingType> >& contactCouplings) {
+
+ contactCoupling_.resize(contactCouplings.size());
+ for (size_t i=0; i<contactCouplings.size(); i++)
+ contactCoupling_[i] = contactCouplings[i];
+ }
+
+ /** \brief Get the transposed of the mortar transformation matrix B^T.*/
+ const MatrixType& getTransformationMatrix() const {return BT_;}
+
+ /** \brief Get the grids. */
+ const std::vector<const GridType*> getGrids() const { return grids_; }
+
+ /** \brief Set the grids. */
+ void setGrids(const std::vector<const GridType*>& grids) {grids_ = grids;}
+
+ /** \brief Set grids. */
+ void setCoupling(const CouplingPair<GridType, GridType, field_type>& coupling, size_t i) {
+ coupling_[i] = coupling;
+ }
+
+ /** \brief Get reference to i'th coupling. */
+ const auto& getCoupling(size_t i) const {return coupling_[i];}
+
+protected:
+ /** \brief Compute the transposed mortar transformation matrix. */
+ void computeTransformationMatrix();
+
+ /** \brief Setup the Householder reflections. */
+ void assembleHouseholderReflections();
+
+public:
+ /** \brief Bitvector for all bodies whose flags are set if a dof has an obstacle.*/
+ Dune::BitSetVector<dim> totalHasObstacle_;
+
+ /** \brief Obstacles for all bodies on the leafview.*/
+ std::vector<BoxConstraint<field_type,dim> > totalObstacles_;
+
+ /** \brief The mortar couplings between grid pairs */
+ std::vector<std::shared_ptr<CouplingType> > contactCoupling_;
+
+ /** \brief The coupling pairs. */
+ std::vector<CouplingPair<GridType,GridType,field_type> > coupling_;
+
+ /** \brief Vector containing the involved grids. */
+ std::vector<const GridType*> grids_;
+
+ /** \brief The transposed of the mortar transformation matrix. */
+ MatrixType BT_;
+
+ /** \brief Dismiss all faces that are not at least covered by the grid-glue projection for this
+ * much percentage ranging between one - for total coverage and zero for taking all faces.
+ */
+ field_type coveredArea_;
+};
+
+} /* namespace Contact */
+} /* namespace Dune */
+
+#include "nbodyassembler.cc"
+
+#endif
diff --git a/dune/tectonic/spatial-solving/preconditioners/nbodycontacttransfer.cc b/dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc
similarity index 100%
rename from dune/tectonic/spatial-solving/preconditioners/nbodycontacttransfer.cc
rename to dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc
diff --git a/dune/tectonic/spatial-solving/preconditioners/nbodycontacttransfer.hh b/dune/tectonic/spatial-solving/contact/nbodycontacttransfer.hh
similarity index 100%
rename from dune/tectonic/spatial-solving/preconditioners/nbodycontacttransfer.hh
rename to dune/tectonic/spatial-solving/contact/nbodycontacttransfer.hh
diff --git a/dune/tectonic/spatial-solving/preconditioners/CMakeLists.txt b/dune/tectonic/spatial-solving/preconditioners/CMakeLists.txt
index bcff7b9ff1be92af3c892c7ac15e61a368a73c23..b32277edecfc293e5d53c4f9f0eccf7ee0f3ffc6 100644
--- a/dune/tectonic/spatial-solving/preconditioners/CMakeLists.txt
+++ b/dune/tectonic/spatial-solving/preconditioners/CMakeLists.txt
@@ -3,8 +3,6 @@ add_custom_target(tectonic_dune_spatial-solving_preconditioners SOURCES
levelpatchpreconditioner.hh
patchproblem.hh
multilevelpatchpreconditioner.hh
- nbodycontacttransfer.hh
- nbodycontacttransfer.cc
supportpatchfactory.hh
)
@@ -14,6 +12,5 @@ install(FILES
levelpatchpreconditioner.hh
patchproblem.hh
multilevelpatchpreconditioner.hh
- nbodycontacttransfer.hh
supportpatchfactory.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/spatial-solving/tnnmg/functional.hh b/dune/tectonic/spatial-solving/tnnmg/functional.hh
index eff8e22d34a8572122b67b503ed573463117d2d5..3e9eab1a55bce71384bd547761cc2501cd30a74c 100755
--- a/dune/tectonic/spatial-solving/tnnmg/functional.hh
+++ b/dune/tectonic/spatial-solving/tnnmg/functional.hh
@@ -422,10 +422,28 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
maxEigenvalues_(this->linearPart().size())
{
for (size_t i=0; i<this->quadraticPart().N(); ++i) {
+ const auto& Aii = this->quadraticPart()[i][i];
+
+ // due to numerical imprecision, Aii might not be symmetric leading to complex eigenvalues
+ // consider Toeplitz decomposition of Aii and take only symmetric part
+
+ auto symBlock = Aii;
+
+ /*for (size_t j=0; j<symBlock.N(); j++) {
+ for (size_t k=0; k<symBlock.M(); k++) {
+ if (symBlock[j][k]/symBlock[j][j] < 1e-14)
+ symBlock[j][k] = 0;
+ }
+ }*/
+
+ try {
typename Vector::block_type eigenvalues;
- Dune::FMatrixHelp::eigenValues(this->quadraticPart()[i][i], eigenvalues);
+ Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
maxEigenvalues_[i] =
*std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+ } catch (Dune::Exception &e) {
+ std::cout << "complex eig in dof: " << i << std::endl;
+ }
}
}
diff --git a/src/foam/foam-2D.cfg b/src/foam/foam-2D.cfg
index 6efcb38b3aa3637d41b65307d5c97f2cb2136428..c22b8a7224393ab58613ba0f4a075e142f412bfc 100644
--- a/src/foam/foam-2D.cfg
+++ b/src/foam/foam-2D.cfg
@@ -1,6 +1,9 @@
# -*- mode:conf -*-
-[boundary.friction]
-smallestDiameter = 0.005 # 2e-3 [m]
+[body0]
+smallestDiameter = 0.01 # 2e-3 [m]
+
+[body1]
+smallestDiameter = 0.01 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index fe95424f3676d3797d2a0b3549a57870d8216edc..ce0f8d3a9ab4f1c188ea71afbb1b4dc9853d3d48 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -330,18 +330,18 @@ int main(int argc, char *argv[]) {
BoundaryNodes dirichletNodes;
contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+ for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
}
- }*/
+ }
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
contactNetwork.frictionNodes(frictionNodes);
- /*for (size_t i=0; i<frictionNodes.size(); i++) {
+ for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
- }*/
+ }
Updaters current(
initRateUpdater(
diff --git a/src/foam/foam.cfg b/src/foam/foam.cfg
index 7c8016105e83a5a03566458f4be8128cd5a12a4a..ee78e2fc57ca31b4e0f2eca7f668fb5098664e29 100644
--- a/src/foam/foam.cfg
+++ b/src/foam/foam.cfg
@@ -2,34 +2,34 @@
gravity = 9.81 # [m/s^2]
[body0]
-length = 0.4 # [m]
+length = 0.04 # [m]
height = 0.04 # [m]
depth = 0.04 # [m]
-bulkModulus = 2190 # [Pa]
-poissonRatio = 0.11 # [1]
+bulkModulus = 0.5e5 # [Pa] #2190
+poissonRatio = 0.30 # [1] #0.11
[body0.elastic]
-density = 750 # [kg/m^3]
-shearViscosity = 0 # [Pas]
-bulkViscosity = 0 # [Pas]
+density = 900 # [kg/m^3] #750
+shearViscosity = 1e3 # [Pas]
+bulkViscosity = 1e3 # [Pas]
[body0.viscoelastic]
-density = 750 # [kg/m^3]
-shearViscosity = 0 # [Pas]
-bulkViscosity = 0 # [Pas]
+density = 1000 # [kg/m^3]
+shearViscosity = 1e4 # [Pas]
+bulkViscosity = 1e4 # [Pas]
[body1]
length = 0.04 # [m]
height = 0.04 # [m]
depth = 0.04 # [m]
-bulkModulus = 2190 # [Pa]
-poissonRatio = 0.11 # [1]
+bulkModulus = 0.5e5 # [Pa]
+poissonRatio = 0.30 # [1]
[body1.elastic]
-density = 750 # [kg/m^3]
-shearViscosity = 0 # [Pas]
-bulkViscosity = 0 # [Pas]
+density = 900 # [kg/m^3]
+shearViscosity = 1e3 # [Pas]
+bulkViscosity = 1e3 # [Pas]
[body1.viscoelastic]
-density = 750 # [kg/m^3]
-shearViscosity = 0 # [Pas]
-bulkViscosity = 0 # [Pas]
+density = 1000 # [kg/m^3]
+shearViscosity = 1e4 # [Pas]
+bulkViscosity = 1e4 # [Pas]
[boundary.friction]
C = 10 # [Pa]
@@ -40,10 +40,10 @@ initialAlpha = 0 # [ ]
stateModel = AgeingLaw
frictionModel = Truncated #Regularised
[boundary.friction.weakening]
-a = 0.025 # [ ]
-b = 0.005 # [ ]
+a = 0.002 # [ ]
+b = 0.017 # [ ]
[boundary.friction.strengthening]
-a = 0.025 # [ ]
+a = 0.020 # [ ]
b = 0.005 # [ ]
[boundary.neumann]