diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7c012be4a60057d74a9d619b34f13b8f3be5ce80..0cfcc7081672871f755ca84b6d7f8d0e39f2e082 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory("tests")
+add_subdirectory("spatial-solving")
set(SW_SOURCE_FILES
assemblers.cc
diff --git a/src/data-structures/levelcontactnetwork.hh b/src/data-structures/levelcontactnetwork.hh
index 97792f2ee6dcd3ff2483e81e1eba0d053a1b3322..68bb716aaac8391fdb2c17798a512c66c91549d7 100644
--- a/src/data-structures/levelcontactnetwork.hh
+++ b/src/data-structures/levelcontactnetwork.hh
@@ -25,6 +25,7 @@
template <class GridType, int dimension> class LevelContactNetwork {
private:
+ //static const int dimension = typename GridType::dimension;
using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
public:
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 1eceb99a8ebdbbd5424d25146e24ce8935dd882d..bab650b970bdd1dffb466e5d648fe0cd3081a4f7 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -16,7 +16,6 @@
#include "levelcontactnetwork.hh"
#include "matrices.hh"
#include "../spatial-solving/solverfactory.hh"
-#include "../spatial-solving/solverfactory_old.hh"
#include "../spatial-solving/tnnmg/functional.hh"
#include "../spatial-solving/tnnmg/zerononlinearity.hh"
#include "../utils/debugutils.hh"
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 5790cc65dfe920c95bda26ee52a3b91ce06d7ebb..e20cc637dd7c631136ba0712c78c3fb8ddc58ead 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -42,8 +42,6 @@
#include <dune/contact/common/deformedcontinuacomplex.hh>
#include <dune/contact/common/couplingpair.hh>
-#include <dune/contact/assemblers/nbodyassembler.hh>
-//#include <dune/contact/assemblers/dualmortarcoupling.hh>
#include <dune/contact/projections/normalprojection.hh>
#include <dune/tectonic/geocoordinate.hh>
@@ -73,6 +71,7 @@
#include "spatial-solving/tnnmg/functional.hh"
#include "spatial-solving/tnnmg/localbisectionsolver.hh"
+#include "spatial-solving/contact/nbodyassembler.hh"
#include "spatial-solving/solverfactory.hh"
#include "time-stepping/adaptivetimestepper.hh"
diff --git a/src/spatial-solving/CMakeLists.txt b/src/spatial-solving/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1970afe64c1da2c2152ba9de21234515d17472fc
--- /dev/null
+++ b/src/spatial-solving/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory("tnnmg")
+add_subdirectory("contact")
+add_subdirectory("preconditioners")
diff --git a/src/spatial-solving/contact/CMakeLists.txt b/src/spatial-solving/contact/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..05ec76c58ad389b8aac6c0368f6f551825f2ccb3
--- /dev/null
+++ b/src/spatial-solving/contact/CMakeLists.txt
@@ -0,0 +1,6 @@
+add_custom_target(dune-tectonic_spatial-solving_contact_src SOURCES
+ dualmortarcoupling.cc
+ dualmortarcoupling.hh
+ nbodyassembler.cc
+ nbodyassembler.hh
+)
diff --git a/src/spatial-solving/contact/dualmortarcoupling.cc b/src/spatial-solving/contact/dualmortarcoupling.cc
new file mode 100644
index 0000000000000000000000000000000000000000..96ca666c521a49612d334915f26adf359af422fa
--- /dev/null
+++ b/src/spatial-solving/contact/dualmortarcoupling.cc
@@ -0,0 +1,333 @@
+// -*- 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>
+
+
+template <class field_type, class GridView0, class GridView1>
+void DualMortarCoupling<field_type, GridView0, GridView1>::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 = gridView0_.indexSet();
+
+ // 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 GridView0, class GridView1>
+void DualMortarCoupling<field_type, GridView0, GridView1>::setup()
+{
+ 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);
+ };
+
+ 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
+ BoundaryPatch0 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();
+}
diff --git a/src/spatial-solving/contact/dualmortarcoupling.hh b/src/spatial-solving/contact/dualmortarcoupling.hh
new file mode 100644
index 0000000000000000000000000000000000000000..90bc49dc759418d2f596a17ec75bc05a8df6b7cf
--- /dev/null
+++ b/src/spatial-solving/contact/dualmortarcoupling.hh
@@ -0,0 +1,199 @@
+// -*- 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>
+
+
+/** \brief Extension of Dune::Contact::DualMortarCoupling that works on
+ * any GridView (not just LeafGridView) and allows crosspoints.
+ *
+ * Assembles the transfer operator for two-body contact
+ */
+template<class field_type, class GridView0, class GridView1=GridView0>
+class DualMortarCoupling {
+
+ // /////////////////////
+ // private types
+ // /////////////////////
+
+ enum {dim = GridView0::dimension};
+
+ enum {dimworld = GridView0::dimensionworld};
+
+ protected:
+ using MatrixBlock = Dune::FieldMatrix<field_type, dim, dim>;
+ using MatrixType = Dune::BCRSMatrix<MatrixBlock>;
+ using ObstacleVector = Dune::BlockVector<Dune::FieldVector<field_type, 1> >;
+
+ using ctype = typename GridView0::ctype;
+
+ using GridType0 = typename GridView0::Grid;
+ using GridType1 = typename GridView1::Grid;
+
+ using BoundaryPatch0 = BoundaryPatch<GridView0>;
+ using BoundaryPatch1 = BoundaryPatch<GridView1>;
+
+ using LevelBoundaryPatch0 = BoundaryPatch<typename GridType0::LevelGridView>;
+ using LevelBoundaryPatch1 = BoundaryPatch<typename GridType1::LevelGridView>;
+
+ using Extractor0 = Dune::GridGlue::Codim1Extractor<GridView0>;
+ using Extractor1 = Dune::GridGlue::Codim1Extractor<GridView1>;
+
+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(GridView0&& gridView0, GridView0&& gridView1,
+ field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
+ : overlap_(overlap), gridGlueBackend_(nullptr), coveredArea_(coveredArea)
+ {
+ setGrids(gridView0, gridView1);
+ }
+
+ 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 nonmortar boundary. */
+ const BoundaryPatch0& nonmortarBoundary() const {
+ return nonmortarBoundary_;
+ }
+
+ /** \brief Return the mortar boundary. */
+ const BoundaryPatch1& mortarBoundary() const {
+ return mortarBoundary_;
+ }
+
+ /** \brief Set the non-mortar and mortar grid. */
+ void setGrids(GridView0&& gridView0, GridView1&& gridView1) {
+ gridView0_ = std::forward<GridView0>(gridView0);
+ gridView1_ = std::forward<GridView1>(gridView1);
+
+ grid0_ = &gridView0_.grid();
+ grid1_ = &gridView1_.grid();
+ }
+
+ /** \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
+ // /////////////////////////////////////////////
+private:
+ /** \brief The two grids involved in the two-body contact problem. */
+ const GridView0& gridView0_;
+ const GridView1& gridView1_;
+
+ const GridType0* grid0_;
+ const GridType1* grid1_;
+
+ /** \brief For each dof a bit specifying whether the dof carries an obstacle or not. */
+ BoundaryPatch0 nonmortarBoundary_;
+
+ /** \brief The mortar boundary. */
+ BoundaryPatch1 mortarBoundary_;
+
+ /** \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_;
+};
+
+
+#include "dualmortarcoupling.cc"
+
+#endif
diff --git a/src/spatial-solving/contact/nbodyassembler.cc b/src/spatial-solving/contact/nbodyassembler.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d8ee25a45d13da8fdaa6f5810937c09631388b64
--- /dev/null
+++ b/src/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/src/spatial-solving/contact/nbodyassembler.hh b/src/spatial-solving/contact/nbodyassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ae1db365287d1355a56753b509f3d6bc4e717395
--- /dev/null
+++ b/src/spatial-solving/contact/nbodyassembler.hh
@@ -0,0 +1,200 @@
+// -*- 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/contactassembler.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/contactprojection.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
+#include "dualmortarcoupling.hh"
+
+#ifdef HAVE_DUNE_GRID_GLUE
+#include <dune/grid-glue/merging/merger.hh>
+#endif
+
+/** \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 Dune::Contact::ContactAssembler<GridType::dimension, typename VectorType::field_type>
+{
+protected:
+ enum {dim = GridType::dimension};
+
+ using field_type = typename Dune::PromotionTraits<typename VectorType::field_type,
+ typename GridType::ctype>::PromotedType;
+
+ using CouplingType = DualMortarCoupling<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_;
+};
+
+#include "nbodyassembler.cc"
+
+#endif
diff --git a/src/spatial-solving/preconditioners/CMakeLists.txt b/src/spatial-solving/preconditioners/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..abb4ff76bea86b3ff68d31eb1b25e92947d747e7
--- /dev/null
+++ b/src/spatial-solving/preconditioners/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_custom_target(dune-tectonic_spatial-solving_preconditioners_src SOURCES
+ hierarchicleveliterator.hh
+ levelglobalpreconditioner.hh
+ levelpatchpreconditioner.hh
+ multilevelpatchpreconditioner.hh
+ supportpatchfactory.hh
+)
diff --git a/src/spatial-solving/preconditioners/hierarchicleveliterator.hh b/src/spatial-solving/preconditioners/hierarchicleveliterator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c8c7cf4cefecc371f30252c2c12c7c22b42159e7
--- /dev/null
+++ b/src/spatial-solving/preconditioners/hierarchicleveliterator.hh
@@ -0,0 +1,155 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef HIERARCHIC_LEVEL_ITERATOR_HH
+#define HIERARCHIC_LEVEL_ITERATOR_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/iteratorfacades.hh>
+
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+
+/** \brief Hierarchic leaf iterator.
+ *
+ * This iterator loops over all children of a given coarse element and only returns the ones that are leaf.
+ * If the starting element is leaf itself, then the returned iterator returns the element itself.
+ * This class also provides a geometry, mapping local coordinates of the children to local coordinates
+ * in the coarse element.
+*/
+template <class GridImp>
+class HierarchicLevelIterator :
+ public Dune::ForwardIteratorFacade<HierarchicLevelIterator<GridImp>, const typename GridImp::template Codim<0>::Entity>
+{
+ typedef typename GridImp::template Codim<0>::Entity Element;
+ typedef typename GridImp::HierarchicIterator HierarchicIterator;
+ enum {dim = GridImp::dimension};
+ enum {dimworld = GridImp::dimensionworld};
+public:
+ typedef Dune::CachedMultiLinearGeometry<typename GridImp::ctype, dim, dimworld> LocalGeometry;
+
+ enum PositionFlag {begin, end};
+
+ HierarchicLevelIterator(const GridImp& grid, const Element& element,
+ PositionFlag flag, int maxlevel, bool nested = true)
+ : element_(element), maxlevel_(maxlevel), hIt_(element.hend(maxlevel_)),
+ flag_(flag), nested_(nested)
+ {
+
+ // if the element itself is leaf, then we don't have to iterate over the children
+ if (flag==begin && !(element_.level()==maxlevel_)) {
+ hIt_ = element_.hbegin(maxlevel_);
+
+ //NOTE This class by now assumes that possible non-nestedness of the grid levels only arises
+ // due to boundary parametrisation
+ if (!nested_) {
+ // Check if the element is a boundary element, and set the nested flag correspondingly
+ // If the element is not at the boundary, then we can neglect the nested flag
+ typename GridImp::LevelIntersectionIterator lIt = grid.levelGridView(element.level()).ibegin(element);
+ typename GridImp::LevelIntersectionIterator lEndIt = grid.levelGridView(element.level()).ibegin(element);
+ nested_ = true;
+ for (; lIt != lEndIt; ++lIt)
+ if (lIt->boundary()) {
+ nested_ = false;
+ break;
+ }
+ }
+
+ if (!(hIt_->level()==maxlevel_))
+ increment();
+ }
+ }
+
+ //! Copy constructor
+ HierarchicLevelIterator(const HierarchicLevelIterator<GridImp>& other) :
+ element_(other.element_),maxlevel_(other.maxlevel_), hIt_(other.hIt_),
+ flag_(other.flag_), nested_(other.nested_)
+ {}
+
+ //! Equality
+ bool equals(const HierarchicLevelIterator<GridImp>& other) const {
+ return (element_ == other.element_) &&
+ ((flag_==other.flag_ && flag_==end) || (hIt_ == other.hIt_ && flag_==begin && other.flag_==flag_));
+ }
+
+ //! Prefix increment
+ void increment() {
+
+ HierarchicIterator hEndIt = element_.hend(maxlevel_);
+
+ if (hIt_ == hEndIt) {
+ flag_ = end;
+ return;
+ }
+
+ // Increment until we hit a leaf child element
+ do {
+ ++hIt_;
+ } while(hIt_ != hEndIt && (!(hIt_->level()==maxlevel_)));
+
+ if (hIt_ == hEndIt)
+ flag_ = end;
+ }
+
+ //! Dereferencing
+ const Element& dereference() const {
+ if (flag_ == end)
+ DUNE_THROW(Dune::Exception,"HierarchicLevelIterator: Cannot dereference end iterator!");
+ return (element_.level()==maxlevel_) ? element_ : *hIt_;
+ }
+
+ //! Compute the local geometry mapping from the leaf child to the starting element
+ LocalGeometry geometryInAncestor() {
+
+ //NOTE: We assume here that the type of the child and the ancestors are the same!
+ const Dune::ReferenceElement<typename GridImp::ctype, dim>& ref = Dune::ReferenceElements<typename GridImp::ctype,dim>::general(element_.type());
+
+ // store local coordinates of the leaf child element within the coarse starting element
+ std::vector<Dune::FieldVector<typename GridImp::ctype,dim> > corners(ref.size(dim));
+ for (int i=0; i<corners.size(); i++)
+ corners[i] = ref.position(i,dim);
+
+ // If the element is leaf, then return an Identity mapping
+ if (element_.level()==maxlevel_)
+ return LocalGeometry(ref,corners);
+
+ if (nested_) {
+ const typename Element::Geometry leafGeom = hIt_->geometry();
+ const typename Element::Geometry coarseGeom = element_.geometry();
+
+ for (int i=0; i<corners.size();i++)
+ corners[i] = coarseGeom.local(leafGeom.corner(i));
+ return LocalGeometry(ref,corners);
+ }
+
+ Element father(*hIt_);
+ while (father != element_) {
+
+ const typename Element::LocalGeometry fatherGeom = hIt_->geometryInFather();
+ father = father->father();
+
+ for (int i=0; i<corners.size(); i++)
+ corners[i] = fatherGeom.global(corners[i]);
+ }
+
+ return LocalGeometry(ref,corners);
+ }
+
+private:
+
+ //! The starting element
+ const Element element_;
+
+ //! The grids maxlevel
+ int maxlevel_;
+
+ //! The actual hierarchic iterator
+ HierarchicIterator hIt_;
+
+ //! Position flag
+ PositionFlag flag_;
+
+ //! Flag that is set true if the grid levels are conforming (i.e. no parametrised boundaries)
+ bool nested_;
+
+};
+#endif
diff --git a/src/spatial-solving/tnnmg/levelglobalpreconditioner.hh b/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
similarity index 98%
rename from src/spatial-solving/tnnmg/levelglobalpreconditioner.hh
rename to src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
index 2b7b4e2bb80f6ad5ec302097654dbb5bc320deb4..ac9fd3fce63579595d9e64d29c2e9e69c1d6b63d 100644
--- a/src/spatial-solving/tnnmg/levelglobalpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
@@ -7,6 +7,7 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+//#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
#include <dune/fufem/boundarypatch.hh>
@@ -115,7 +116,7 @@ class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorT
this->x_->resize(matrix_.M());
*(this->x_) = 0;
- #if HAVE_UMFPACK
+ #if HAVE_SUPERLU
Dune::InverseOperatorResult res;
VectorType rhsCopy(rhs_);
@@ -123,7 +124,7 @@ class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorT
solver.apply(*(this->x_), rhsCopy, res);
#else
- #error No UMFPack!
+ #error No SuperLU!
#endif
//std::cout << "LevelGlobalPreconditioner::iterate() Solving global problem took: " << timer.elapsed() << " seconds" << std::endl;
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d905e4a3b40f47cba638278ee36706ea88695ba5
--- /dev/null
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -0,0 +1,210 @@
+#ifndef LEVEL_PATCH_PRECONDITIONER_HH
+#define LEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include "../../data-structures/levelcontactnetwork.hh"
+
+#include "supportpatchfactory.hh"
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+
+template <class PatchSolver, class MatrixType, class VectorType>
+class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+ enum Mode {additive, multiplicative};
+
+private:
+ using Base = LinearIterationStep<MatrixType, VectorType>;
+
+ using PatchFactory = SupportPatchFactory<GridView, GridView>;
+ using Patch = typename SupportPatchFactory::Patch;
+
+ typedef typename BasisType::GridView GridView;
+ typedef typename GridView::Grid GridType;
+
+ static const int dim = GridType::dimension;
+
+ const Mode mode_;
+
+ const LevelContactNetwork<GridView, dim>& levelContactNetwork_;
+
+ std::vector<Patch> patches_;
+ PatchSolver patchSolver_;
+
+ const GridType& grid_;
+ const int level_;
+ const BasisType basis_;
+
+ size_t patchDepth_;
+
+ MatrixType matrix_;
+ std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
+
+public:
+
+ // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
+ LevelPatchPreconditioner(const LevelContactNetwork<GridView, dim>& levelContactNetwork,
+ const BasisType& patchLevelBasis,
+ const LocalAssembler& localAssembler,
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ const Mode mode = LevelPatchPreconditioner::Mode::additive) :
+ levelInterfaceNetwork_(levelInterfaceNetwork),
+ patchLevelBasis_(patchLevelBasis),
+ localAssembler_(localAssembler),
+ localInterfaceAssemblers_(localInterfaceAssemblers),
+ mode_(mode),
+ grid_(levelInterfaceNetwork_.grid()),
+ level_(levelInterfaceNetwork_.level()),
+ basis_(levelInterfaceNetwork_)
+ {
+ setPatchDepth();
+ this->verbosity_ = QUIET;
+ }
+
+ // build support patches
+ void build() {
+
+
+ localProblems_.resize(vertexInElements.size());
+
+ // init local fine level corrections
+ if (this->verbosity_ == FULL) {
+ std::cout << std::endl;
+ std::cout << "---------------------------------------------" << std::endl;
+ std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
+
+ Dune::Timer timer;
+ timer.reset();
+ timer.start();
+ }
+
+ const int patchLevel = patchLevelBasis_.faultNetwork().level();
+ for (size_t i=0; i<vertexInElements.size(); i++) {
+ if (this->verbosity_ != QUIET) {
+ std::cout << i << std::endl;
+ std::cout << "---------------" << std::endl;
+ }
+
+ Patch patchDofs;
+ SupportPatchFactory<decltype(coarseGridView), decltype(fineGridView)> patchFactory(coarseGridView, fineGridView, vertexInElementsBody);
+ patchFactory.build(seedObject, patchDofs, patchDepth_);
+
+ if ((i+1) % 10 == 0 && this->verbosity_ == FULL) {
+ std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
+ }
+ }
+
+ if (this->verbosity_ == FULL) {
+ std::cout << std::endl;
+ std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
+ std::cout << "---------------------------------------------" << std::endl;
+ timer.stop();
+ }
+
+ if (this->verbosity_ != REDUCED) {
+ std::cout << "Level: " << level_ << " LevelPatchPreconditioner::build() - SUCCESS" << std::endl;
+ }
+ }
+
+ void setPatchDepth(const size_t patchDepth = 0) {
+ patchDepth_ = patchDepth;
+ }
+
+ void setVerbosity(VerbosityMode verbosity) {
+ this->verbosity_ = verbosity;
+ }
+
+ virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
+
+ /* VectorType rhs;
+ rhs.resize(matrix_.N());
+ rhs = 0;
+
+ //print(localToGlobal, "localToGlobal: ");
+ //print(boundaryDofs, "boundaryDofs: ");
+
+ localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);*/
+
+ this->x_ = &x;
+ this->rhs_ = &rhs;
+ this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+
+ patchSolver_.setProblem(mat, x, rhs);
+ }
+
+ virtual void iterate() {
+ if (mode_ == additive)
+ iterateAdd();
+ else
+ iterateMult();
+ }
+
+ void iterateAdd() {
+ *(this->x_) = 0;
+
+ VectorType x = *(this->x_);
+ for (size_t i=0; i<patches_.size(); i++) {
+ patchSolver_.setIgnore(patches_[i]);
+ patchSolver_.solve();
+
+ *(this->x_) += x;
+ }
+ }
+
+ void iterateMult() {
+ *(this->x_) = 0;
+
+ VectorType x = *(this->x_);
+ for (size_t i=0; i<patches_.size(); i++) {
+ VectorType updatedRhs(*(this->rhs_));
+ this->mat_->mmv(*(this->x_), updatedRhs);
+
+ //step(i);
+ //print(updatedRhs, "localRhs: ");
+ //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));
+
+ patchSolver_.setProblem(*this->mat_, *(this->x_), updatedRhs);
+ patchSolver_.setIgnore(patches_[i]);
+ patchSolver_.solve();
+
+ //print(it, "it: ");
+ //print(x, "x: ");
+
+ //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));
+
+ /*if (i==5) {
+ writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
+ }*/
+
+ *(this->x_) += x;
+ }
+ }
+
+ const GridView& gridView() const {
+ return basis_.getGridView();
+ }
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
+ return levelInterfaceNetwork_;
+ }
+
+ size_t size() const {
+ return localProblems_.size();
+ }
+
+};
+
+#endif
+
diff --git a/src/spatial-solving/tnnmg/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
similarity index 100%
rename from src/spatial-solving/tnnmg/multilevelpatchpreconditioner.hh
rename to src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f518ecc0de92b9fcdbaab92bdfca58b36f26072a
--- /dev/null
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -0,0 +1,514 @@
+#ifndef SUPPORT_PATCH_FACTORY_HH
+#define SUPPORT_PATCH_FACTORY_HH
+
+#include <set>
+#include <queue>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+//#include "../../data-structures/levelcontactnetwork.hh"
+#include "hierarchicleveliterator.hh"
+
+template <class Element>
+class BodyElement {
+private:
+ size_t bodyID_;
+ Element element_;
+
+public:
+ BodyElement() {}
+ BodyElement(size_t bodyID, const Element& element) {
+ set(bodyID, element);
+ }
+
+ void set(size_t bodyID, const Element& element) {
+ bodyID_ = bodyID;
+ element_ = element;
+ }
+
+ const auto& bodyID() const {
+ return bodyID_;
+ }
+ const auto& element() const {
+ return element_;
+ }
+};
+
+
+template <class LevelContactNetwork>
+class BodyHelper {
+private:
+ using GridType = typename LevelContactNetwork::DeformedGridType;
+ static const int dim = GridType::dimension;
+
+ using Element = typename GridType::template Codim<0>::Entity;
+ //using Vertex = typename GridType::template Codim<dim>::Entity;
+
+ using LevelGridView = typename LevelContactNetwork::DeformedLevelGridView;
+ using LeafGridView = typename LevelContactNetwork::DeformedLeafGridView;
+
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<GridType>;
+ using LevelElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, Dune::MCMGElementLayout >;
+ using LevelFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, mcmgLayout(Codim<1>()) >;
+
+ using LeafElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, Dune::MCMGElementLayout >;
+ using LeafFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, mcmgLayout(Codim<1>()) >;
+
+ bool levelMode_;
+
+ std::unique_ptr<LevelGridView> levelGridView_;
+ LevelElementMapper levelElementMapper_;
+ LevelFaceMapper levelFaceMapper_;
+
+ std::unique_ptr<LeafGridView> leafGridView_;
+ LeafElementMapper leafElementMapper_;
+ LeafFaceMapper leafFaceMapper_;
+
+ std::vector<std::vector<Element>>& vertexInElements_;
+
+ //std::vector<size_t> globalToLocal_;
+
+ template <class GV>
+ void init(const GV& gv) {
+ // init vertexInElements
+ vertexInElements_.resize(gv.size(dim));
+ for (size_t i=0; i<vertexInElements_.size(); i++) {
+ vertexInElements_[i].resize(0);
+ }
+
+ const auto& indexSet = gv.indexSet();
+
+ for (const auto& elem : elements(gv)) {
+ const auto& lfe = cache_.get(elem.type());
+
+ for (size_t j=0; j<lfe.localBasis().size(); ++j) {
+ auto localIdx = lfe.localCoefficients().localKey(j).subEntity();
+ int dofIndex = indexSet.subIndex(elem, localIdx, dim);
+
+ vertexInElements_[dofIndex].push_back(elem);
+ }
+ }
+ }
+public:
+
+ BodyHelper(const LevelGridView& gridView) :
+ levelMode_(true),
+ levelGridView_(std::make_unique<LevelGridView>(gridView)),
+ levelElementMapper_(gridView),
+ leafFaceMapper_(gridView) {
+
+ init(gridView);
+ }
+
+ BodyHelper(const LeafGridView& gridView) :
+ levelMode_(false),
+ leafGridView_(std::make_unique<LeafGridView>(gridView)),
+ leafElementMapper_(gridView),
+ leafFaceMapper_(gridView) {
+
+ init(gridView);
+ }
+
+ const auto& gridView() const {
+ return (levelMode_ ? levelGridView_ : leafGridView_);
+ }
+
+ const auto& elementMapper() const {
+ return (levelMode_ ? levelElementMapper_ : leafElementMapper_);
+ }
+
+ const auto& faceMapper() const {
+ return (levelMode_ ? levelFaceMapper_ : leafFaceMapper_);
+ }
+
+ const auto& vertexInElements() const {
+ return vertexInElements_;
+ }
+
+ const auto& cache() const {
+ return cache_;
+ }
+};
+
+
+/**
+ * constructs BitSetVector whose entries are
+ * true, if vertex is outside of patch or part of patch boundary
+ * false, if vertex is in interior of patch
+ */
+template <class LevelContactNetwork>
+class SupportPatchFactory
+{
+public:
+ using Patch = Dune::BitSetVector<1>;
+
+private:
+ static const int dim = LevelContactNetwork::dimension; //TODO
+ using ctype = typename LevelContactNetwork::ctype;
+
+ using BodyHelper = BodyHelper<LevelContactNetwork>;
+ //using Vertex = typename GridType::template Codim<dim>::Entity;
+
+ using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity;
+ using BodyElement = BodyElement<Element>;
+ using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
+
+ const size_t nBodies_;
+
+ const LevelContactNetwork& coarseContactNetwork_;
+ const LevelContactNetwork& fineContactNetwork_;
+
+ std::vector<std::unique_ptr<BodyHelper>> bodyHelper_;
+
+ std::map<size_t, std::vector<BodyElement>> neighborElements_; // map faceID to neighbor elements
+ std::map<size_t, std::vector<size_t>> neighborElementDofs_; // map faceID to neighbor element dofs
+
+ std::vector<size_t> vertexOffSet_;
+ std::vector<size_t> faceOffSet_;
+ std::vector<size_t> elementOffSet_;
+
+ FiniteElementCache cache_;
+
+public:
+ SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
+ nBodies_(coarseContactNetwork_.nBodies()),
+ coarseContactNetwork_(coarseContactNetwork),
+ fineContactNetwork_(fineContactNetwork),
+ bodyHelper_(nBodies_),
+ vertexOffSet_(nBodies_),
+ faceOffSet_(nBodies_),
+ elementOffSet_(nBodies_){
+
+ assert(nBodies_ == fineContactNetwork_.nBodies());
+
+ size_t vertexOffSet = 0;
+ size_t faceOffSet = 0;
+ size_t elementOffSet = 0;
+
+ for (size_t i=0; i<nBodies_; i++) {
+ const auto& gridView = coarseContactNetwork_.body(i)->gridView();
+ bodyHelper_[i] = std::make_unique<BodyHelper>(gridView);
+
+ vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
+ vertexOffSet = vertexOffSet_[i];
+
+ faceOffSet_[i] = faceOffSet + gridView.size(1);
+ faceOffSet = faceOffSet_[i];
+
+ elementOffSet_[i] = elementOffSet + gridView.size(0);
+ elementOffSet = elementOffSet_[i];
+ }
+
+ // init boundary mapping/transfer
+ const auto& nBodyAssembler = coarseContactNetwork_.nBodyAssembler();
+ const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ const auto& coupling = nBodyAssembler.getCoupling(i);
+ const auto& contactCoupling = contactCouplings[i];
+
+ const size_t nmBody = coupling.gridIdx_[0];
+ const size_t mBody = coupling.gridIdx_[1];
+
+ const auto& nmBoundary = contactCoupling.nonmortarBoundary();
+ const auto& mBoundary = contactCoupling.mortarBoundary();
+
+ const auto& nmIndexSet = nonmortarBoundary.gridView().indexSet();
+ const auto& mIndexSet = mortarBoundary.gridView().indexSet();
+
+ for (const auto& rIs : intersections(*contactCoupling.getGlue())) {
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ int nmFaceIdx = nmIndexSet.subIndex(inside, rIs.indexInInside(), 1);
+ int mFaceIdx = mIndexSet.subIndex(outside, rIs.indexInOutside(), 1);
+
+ BodyElement nmBodyElement(nmBody, inside);
+ BodyElement mBodyElement(mBody, outside);
+
+ neighborElements_[nmFaceIdx].emplace_back(mBodyElement);
+ neighborElements_[mFaceIdx].emplace_back(nmBodyElement);
+
+ std::set<size_t> dofs;
+ computeIntersectionDofs(nmIndexSet, rIs, dofs);
+ neighborElementDofs_[nmFaceIdx].insert(neighborElementDofs_[nmFaceIdx].end(), dofs.begin(), dofs.end());
+
+ computeIntersectionDofs(mIndexSet, rIs, dofs, false);
+ neighborElementDofs_[mFaceIdx].insert(neighborElementDofs_[mFaceIdx].end(), dofs.begin(), dofs.end());
+ }
+ }
+
+ // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
+ // boundary dofs are contained once, interior dofs multiple times
+ // task: remove boundary dofs and redundant multiples of interior dofs
+ for (auto& it : neighborElementDofs_) {
+ auto& dofs = it.second;
+
+ std::sort(dofs.begin(), dofs.end());
+ auto& dof = dofs.begin();
+ while (dof != dofs.end()) {
+ auto next = dof;
+ if (*dof != *(++next))
+ dof = dofs.erase(dof);
+ }
+ dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
+ }
+ }
+
+ void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0) {
+ BodyElement seedElement(bodyID, coarseElement);
+
+ // construct coarse patch
+ std::vector<BodyElement> patchElements;
+ patchElements.emplace_back(seedElement);
+
+ std::set<size_t> visitedElements;
+ visitedElements.insert(elementIndex(seedElement));
+
+ std::set<size_t> patchVertices;
+ patchVertices.insert(vertexIndex(bodyID, coarseElement, localVertex));
+
+ std::queue<BodyElement> patchSeeds;
+ patchSeeds.push(seedElement);
+
+ for (size_t depth=0; depth<=patchDepth; depth++) {
+ std::vector<BodyElement> nextSeeds;
+ std::set<size_t> newPatchVertices;
+
+ while (!patchSeeds.empty()) {
+ const auto& patchSeed = patchSeeds.front();
+ patchSeeds.pop();
+
+ size_t elemIdx = elementIndex(patchSeed);
+
+ if (visitedElements.count(elemIdx))
+ continue;
+
+ visitedElements.insert(elemIdx);
+
+ size_t bodyIdx = patchSeed.bodyID();
+ const auto& elem = patchSeed.element();
+ const auto& gridView = bodyHelper_[bodyIdx]->gridView();
+
+ for (auto& it : intersections(gridView, elem)) {
+ if (isPatchIntersection(gridView.indexSet(), it, patchVertices, newPatchVertices)) {
+ if (it.neighbor()) {
+ BodyElement neighbor(bodyIdx, it.outside());
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ patchElements.emplace_back(neighbor);
+ patchSeeds.push(neighbor);
+ } else {
+ size_t faceIdx = faceIndex(bodyIdx, it);
+
+ if (neighborElements_.count(faceIdx)) {
+ const auto& neighbors = neighborElements_[faceIdx];
+ for (size_t i=0; i<neighbors.size(); i++) {
+ const auto& neighbor = neighbors[i];
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
+ patchElements.emplace_back(neighbor);
+ patchSeeds.push(neighbor);
+ }
+ }
+ }
+ } else {
+ if (it.neighbor()) {
+ BodyElement neighbor(bodyIdx, it.outside());
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ nextSeeds.emplace_back(neighbor);
+ } else {
+ size_t faceIdx = faceIndex(bodyIdx, it);
+
+ if (neighborElements_.count(faceIdx)) {
+ const auto& neighbors = neighborElements_[faceIdx];
+ for (size_t i=0; i<neighbors.size(); i++) {
+ const auto& neighbor = neighbors[i];
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
+ nextSeeds.emplace_back(neighbor);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (size_t i=0; i<nextSeeds.size(); i++) {
+ patchSeeds.push(nextSeeds[i]);
+ }
+ patchVertices.insert(newPatchVertices.begin(), newPatchVertices.end());
+ }
+
+ // construct fine patch
+ for (size_t i=0; i<patchElements.size(); ++i) {
+ addFinePatchElements(patchElements[i], visitedElements, patchDofs);
+ }
+ }
+
+private:
+ template <class IndexSet, class Intersection>
+ bool isPatchIntersection(const IndexSet& indexSet, const Intersection& is, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
+ std::set<size_t> intersectionDofs;
+ computeIntersectionDofs(indexSet, is, intersectionDofs);
+
+ for (auto& dof : intersectionDofs) {
+ if (patchVertices.count(dof)) {
+ newPatchVertices.insert(intersectionDofs.begin(), intersectionDofs.end());
+ return true;
+ }
+ }
+ return false;
+ }
+
+ size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
+ const auto& gridView = bodyHelper_[bodyID]->gridView();
+
+ auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
+ return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
+ }
+
+ template <class Intersection>
+ size_t faceIndex(const size_t bodyID, const Intersection& is) const {
+ const auto& gridView = bodyHelper_[bodyID]->gridView();
+ return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
+ }
+
+ size_t elementIndex(const BodyElement& bodyElem) const {
+ const auto& gridView = bodyHelper_[bodyElem.bodyID()]->gridView();
+ return elementOffSet_[bodyElem.bodyID()] + gridView.indexSet().index(bodyElem.element());
+ }
+
+ template <class Intersection>
+ bool containsInsideSubentity(const Element& elem, const Intersection& intersection, int subEntity, int codim) {
+ return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
+ }
+
+ void addFinePatchElements(const BodyElement& coarseElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
+
+ const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID())->gridView().grid();
+ const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID())->gridView().level();
+
+ HierarchicLevelIteratorType endHierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
+ HierarchicLevelIteratorType hierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
+ for (; hierIt!=endHierIt; ++hierIt) {
+ const Element& fineElem = *hierIt;
+ addLocalDofs(coarseElem, fineElem, coarsePatch, patchDofs);
+ }
+ }
+
+ void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
+ // insert local dofs
+ const auto& lfe = cache_.get(fineElem.type());
+ const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID())->gridView();
+ const auto& indexSet = gridView.indexSet();
+
+ for (size_t i=0; i<lfe.localBasis().size(); ++i) {
+ auto localIdx = lfe.localCoefficients().localKey(i).subEntity();
+ int dofIndex = indexSet.subIndex(fineElem, localIdx, dim);
+ patchDofs[dofIndex][0] = false;
+ }
+
+ const auto coarseLevel = coarseElem.element().level();
+
+ // search for parts of the global and inner boundary
+ for (const auto& is : intersections(gridView, fineElem)) {
+ bool isPatchBoundary = false;
+
+ if (is.neighbor()) {
+ auto outsideFather = is.outside();
+
+ while (outsideFather.level()>coarseLevel) {
+ outsideFather = outsideFather.father();
+ }
+
+ BodyElement outside(coarseElem.bodyID(), outsideFather);
+
+ if (!coarsePatch.count(elementIndex(outside))) {
+ isPatchBoundary = true;
+ }
+ } else {
+ isPatchBoundary = true;
+
+ auto outsideFather = is.outside();
+
+ while (outsideFather.level()>coarseLevel) {
+ outsideFather = outsideFather.father();
+ }
+
+ BodyElement outside(coarseElem.bodyID(), outsideFather);
+
+ size_t faceIdx = faceIndex(coarseElem.bodyID(), is);
+
+ if (neighborElements_.count(faceIdx)) {
+ const auto& outside = neighbors[0];
+
+ if (coarsePatch.count(elementIndex(outside))) {
+ isPatchBoundary = false;
+ }
+ }
+ }
+
+ if (isPatchBoundary) {
+ const auto& localCoefficients = cache_.get(fineElem.type()).localCoefficients();
+
+ for (size_t i=0; i<localCoefficients.size(); i++) {
+ auto entity = localCoefficients.localKey(i).subEntity();
+ auto codim = localCoefficients.localKey(i).codim();
+
+ if (containsInsideSubentity(fineElem, is, entity, codim)) {
+ int dofIndex = indexSet.subIndex(fineElem, entity, dim);
+ patchDofs[dofIndex][0] = true;
+ }
+ }
+ }
+ }
+ }
+
+
+ template <class IndexSet, class Intersection>
+ void computeIntersectionDofs(const IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs, bool inside = true) {
+ intersectionDofs.clear();
+
+ Element* elem;
+ size_t intersectionIndex;
+
+ if (inside) {
+ elem = &intersection.inside();
+ intersectionIndex = intersection.indexInInside();
+ } else {
+ elem = &intersection.outside();
+ intersectionIndex = intersection.indexInOutside();
+ }
+
+ const int dimElement = Element::dimension;
+ const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem->type());
+
+ for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
+ size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
+ size_t globalIdx = idxSet.subIndex(*elem, idxInElement, dimElement);
+
+ intersectionDofs.insert(globalIdx);
+ }
+ }
+};
+#endif
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave b/src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave
new file mode 100644
index 0000000000000000000000000000000000000000..869cbd83d4f0cf01737f5d7d41f7cd53c392da96
--- /dev/null
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave
@@ -0,0 +1,518 @@
+#ifndef SUPPORT_PATCH_FACTORY_HH
+#define SUPPORT_PATCH_FACTORY_HH
+
+#include <set>
+#include <queue>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+//#include "../../data-structures/levelcontactnetwork.hh"
+#include "hierarchicleveliterator.hh"
+
+template <class Element>
+class BodyElement {
+private:
+ size_t bodyID_;
+ Element element_;
+
+public:
+ BodyElement() {}
+ BodyElement(size_t bodyID, const Element& element) {
+ set(bodyID, element);
+ }
+
+ void set(size_t bodyID, const Element& element) {
+ bodyID_ = bodyID;
+ element_ = element;
+ }
+
+ const auto& bodyID() const {
+ return bodyID_;
+ }
+ const auto& element() const {
+ return element_;
+ }
+};
+
+
+template <class LevelContactNetwork>
+class BodyHelper {
+private:
+ using GridType = typename LevelContactNetwork::DeformedGridType;
+ static const int dim = GridType::dimension;
+
+ using Element = typename GridType::template Codim<0>::Entity;
+ //using Vertex = typename GridType::template Codim<dim>::Entity;
+
+ using LevelGridView = typename LevelContactNetwork::DeformedLevelGridView;
+ using LeafGridView = typename LevelContactNetwork::DeformedLeafGridView;
+
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<GridType>;
+ using LevelElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, Dune::MCMGElementLayout >;
+ using LevelFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, mcmgLayout(Codim<1>()) >;
+
+ using LeafElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, Dune::MCMGElementLayout >;
+ using LeafFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, mcmgLayout(Codim<1>()) >;
+
+ bool levelMode_;
+
+ std::unique_ptr<LevelGridView> levelGridView_;
+ LevelElementMapper levelElementMapper_;
+ LevelFaceMapper levelFaceMapper_;
+
+ std::unique_ptr<LeafGridView> leafGridView_;
+ LeafElementMapper leafElementMapper_;
+ LeafFaceMapper leafFaceMapper_;
+
+ std::vector<std::vector<Element>>& vertexInElements_;
+
+ //std::vector<size_t> globalToLocal_;
+
+ template <class GV>
+ void init(const GV& gv) {
+ // init vertexInElements
+ vertexInElements_.resize(gv.size(dim));
+ for (size_t i=0; i<vertexInElements_.size(); i++) {
+ vertexInElements_[i].resize(0);
+ }
+
+ const auto& indexSet = gv.indexSet();
+
+ for (const auto& elem : elements(gv)) {
+ const auto& lfe = cache_.get(elem.type());
+
+ for (size_t j=0; j<lfe.localBasis().size(); ++j) {
+ auto localIdx = lfe.localCoefficients().localKey(j).subEntity();
+ int dofIndex = indexSet.subIndex(elem, localIdx, dim);
+
+ vertexInElements_[dofIndex].push_back(elem);
+ }
+ }
+ }
+public:
+
+ BodyHelper(const LevelGridView& gridView) :
+ levelMode_(true),
+ levelGridView_(std::make_unique<LevelGridView>(gridView)),
+ levelElementMapper_(gridView),
+ leafFaceMapper_(gridView) {
+
+ init(gridView);
+ }
+
+ BodyHelper(const LeafGridView& gridView) :
+ levelMode_(false),
+ leafGridView_(std::make_unique<LeafGridView>(gridView)),
+ leafElementMapper_(gridView),
+ leafFaceMapper_(gridView) {
+
+ init(gridView);
+ }
+
+ const auto& gridView() const {
+ return (levelMode_ ? levelGridView_ : leafGridView_);
+ }
+
+ const auto& elementMapper() const {
+ return (levelMode_ ? levelElementMapper_ : leafElementMapper_);
+ }
+
+ const auto& faceMapper() const {
+ return (levelMode_ ? levelFaceMapper_ : leafFaceMapper_);
+ }
+
+ const auto& vertexInElements() const {
+ return vertexInElements_;
+ }
+
+ const auto& cache() const {
+ return cache_;
+ }
+};
+
+
+/**
+ * constructs BitSetVector whose entries are
+ * true, if vertex is outside of patch or part of patch boundary
+ * false, if vertex is in interior of patch
+ */
+template <class LevelContactNetwork>
+class SupportPatchFactory
+{
+public:
+ using Patch = Dune::BitSetVector<1>;
+
+private:
+ static const int dim = LevelContactNetwork::dimension; //TODO
+ using ctype = typename LevelContactNetwork::ctype;
+
+ using BodyHelper = BodyHelper<LevelContactNetwork>;
+ //using Vertex = typename GridType::template Codim<dim>::Entity;
+
+ using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity;
+ using BodyElement = BodyElement<Element>;
+ using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
+
+ const size_t nBodies_;
+
+ const LevelContactNetwork& coarseContactNetwork_;
+ const LevelContactNetwork& fineContactNetwork_;
+
+ std::vector<std::unique_ptr<BodyHelper>> bodyHelper_;
+
+ std::map<size_t, std::vector<BodyElement>> neighborElements_; // map faceID to neighbor elements
+ std::map<size_t, std::vector<size_t>> neighborElementDofs_; // map faceID to neighbor element dofs
+
+ std::vector<size_t> vertexOffSet_;
+ std::vector<size_t> faceOffSet_;
+ std::vector<size_t> elementOffSet_;
+
+ FiniteElementCache cache_;
+
+public:
+ SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
+ nBodies_(coarseContactNetwork_.nBodies()),
+ coarseContactNetwork_(coarseContactNetwork),
+ fineContactNetwork_(fineContactNetwork),
+ bodyHelper_(nBodies_),
+ vertexOffSet_(nBodies_),
+ faceOffSet_(nBodies_),
+ elementOffSet_(nBodies_){
+
+ assert(nBodies_ == fineContactNetwork_.nBodies());
+
+ size_t vertexOffSet = 0;
+ size_t faceOffSet = 0;
+ size_t elementOffSet = 0;
+
+ for (size_t i=0; i<nBodies_; i++) {
+ const auto& gridView = coarseContactNetwork_.body(i)->gridView();
+ bodyHelper_[i] = std::make_unique<BodyHelper>(gridView);
+
+ vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
+ vertexOffSet = vertexOffSet_[i];
+
+ faceOffSet_[i] = faceOffSet + gridView.size(1);
+ faceOffSet = faceOffSet_[i];
+
+ elementOffSet_[i] = elementOffSet + gridView.size(0);
+ elementOffSet = elementOffSet_[i];
+ }
+
+ // init boundary mapping/transfer
+ const auto& nBodyAssembler = coarseContactNetwork_.nBodyAssembler();
+ const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ const auto& coupling = nBodyAssembler.getCoupling(i);
+ const auto& contactCoupling = contactCouplings[i];
+
+ const size_t nmBody = coupling.gridIdx_[0];
+ const size_t mBody = coupling.gridIdx_[1];
+
+ const auto& nmBoundary = contactCoupling.nonmortarBoundary();
+ const auto& mBoundary = contactCoupling.mortarBoundary();
+
+ const auto& nmIndexSet = nonmortarBoundary.gridView().indexSet();
+ const auto& mIndexSet = mortarBoundary.gridView().indexSet();
+
+ for (const auto& rIs : intersections(*contactCoupling.getGlue())) {
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ int nmFaceIdx = nmIndexSet.subIndex(inside, rIs.indexInInside(), 1);
+ int mFaceIdx = mIndexSet.subIndex(outside, rIs.indexInOutside(), 1);
+
+ BodyElement nmBodyElement(nmBody, inside);
+ BodyElement mBodyElement(mBody, outside);
+
+ neighborElements_[nmFaceIdx].emplace_back(mBodyElement);
+ neighborElements_[mFaceIdx].emplace_back(nmBodyElement);
+
+ std::set<size_t> dofs;
+ computeIntersectionDofs(nmIndexSet, rIs, dofs);
+ neighborElementDofs_[nmFaceIdx].insert(neighborElementDofs_[nmFaceIdx].end(), dofs.begin(), dofs.end());
+
+ computeIntersectionDofs(mIndexSet, rIs, dofs, false);
+ neighborElementDofs_[mFaceIdx].insert(neighborElementDofs_[mFaceIdx].end(), dofs.begin(), dofs.end());
+ }
+ }
+
+ // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
+ // boundary dofs are contained once, interior dofs multiple times
+ // task: remove boundary dofs and redundant multiples of interior dofs
+ for (auto& it : neighborElementDofs_) {
+ auto& dofs = it.second;
+
+ std::sort(dofs.begin(), dofs.end());
+ auto& dof = dofs.begin();
+ while (dof != dofs.end()) {
+ auto next = dof;
+ if (*dof != *(++next))
+ dof = dofs.erase(dof);
+ }
+ dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
+ }
+ }
+
+ void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0) {
+ BodyElement seedElement(bodyID, coarseElement);
+
+ // construct coarse patch
+ std::vector<BodyElement> patchElements;
+ patchElements.emplace_back(seedElement);
+
+ std::set<size_t> visitedElements;
+ visitedElements.insert(elementIndex(seedElement));
+
+ std::set<size_t> patchVertices;
+ patchVertices.insert(vertexIndex(bodyID, coarseElement, localVertex));
+
+ std::queue<BodyElement> patchSeeds;
+ patchSeeds.push(seedElement);
+
+ for (size_t depth=0; depth<=patchDepth; depth++) {
+ std::vector<BodyElement> nextSeeds;
+ std::set<size_t> newPatchVertices;
+
+ while (!patchSeeds.empty()) {
+ const auto& patchSeed = patchSeeds.front();
+ patchSeeds.pop();
+
+ size_t elemIdx = elementIndex(patchSeed);
+
+ if (visitedElements.count(elemIdx))
+ continue;
+
+ visitedElements.insert(elemIdx);
+
+ size_t bodyIdx = patchSeed.bodyID();
+ const auto& elem = patchSeed.element();
+ const auto& gridView = bodyHelper_[bodyIdx]->gridView();
+
+ for (auto& it : intersections(gridView, elem)) {
+ if (isPatchIntersection(gridView.indexSet(), it, patchVertices, newPatchVertices)) {
+ if (it.neighbor()) {
+ BodyElement neighbor(bodyIdx, it.outside());
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ patchElements.emplace_back(neighbor);
+ patchSeeds.push(neighbor);
+ } else {
+ size_t faceIdx = faceIndex(bodyIdx, it);
+
+ if (neighborElements_.count(faceIdx)) {
+ const auto& neighbors = neighborElements_[faceIdx];
+ for (size_t i=0; i<neighbors.size(); i++) {
+ const auto& neighbor = neighbors[i];
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
+ patchElements.emplace_back(neighbor);
+ patchSeeds.push(neighbor);
+ }
+ }
+ }
+ } else {
+ if (it.neighbor()) {
+ BodyElement neighbor(bodyIdx, it.outside());
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ nextSeeds.emplace_back(neighbor);
+ } else {
+ size_t faceIdx = faceIndex(bodyIdx, it);
+
+ if (neighborElements_.count(faceIdx)) {
+ const auto& neighbors = neighborElements_[faceIdx];
+ for (size_t i=0; i<neighbors.size(); i++) {
+ const auto& neighbor = neighbors[i];
+
+ if (visitedElements.count(elementIndex(neighbor)))
+ continue;
+
+ newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
+ nextSeeds.emplace_back(neighbor);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (size_t i=0; i<nextSeeds.size(); i++) {
+ patchSeeds.push(nextSeeds[i]);
+ }
+ patchVertices.insert(newPatchVertices.begin(), newPatchVertices.end());
+ }
+
+ // construct fine patch
+ for (size_t i=0; i<patchElements.size(); ++i) {
+ addFinePatchElements(patchElements[i], visitedElements, patchDofs);
+ }
+ }
+
+private:
+ template <class IndexSet, class Intersection>
+ bool isPatchIntersection(const IndexSet& indexSet, const Intersection& is, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
+ std::set<size_t> intersectionDofs;
+ computeIntersectionDofs(indexSet, is, intersectionDofs);
+
+ for (auto& dof : intersectionDofs) {
+ if (patchVertices.count(dof)) {
+ newPatchVertices.insert(intersectionDofs.begin(), intersectionDofs.end());
+ return true;
+ }
+ }
+ return false;
+ }
+
+ size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
+ const auto& gridView = bodyHelper_[bodyID]->gridView();
+
+ auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
+ return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
+ }
+
+ template <class Intersection>
+ size_t faceIndex(const size_t bodyID, const Intersection& is) const {
+ const auto& gridView = bodyHelper_[bodyID]->gridView();
+ return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
+ }
+
+ size_t elementIndex(const BodyElement& bodyElem) const {
+ const auto& gridView = bodyHelper_[bodyElem.bodyID()]->gridView();
+ return elementOffSet_[bodyElem.bodyID()] + gridView.indexSet().index(bodyElem.element());
+ }
+
+ template <class Intersection>
+ bool containsInsideSubentity(const Element& elem, const Intersection& intersection, int subEntity, int codim) {
+ return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
+ }
+
+ void addFinePatchElements(const BodyElement& coarseElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
+
+ const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID())->gridView().grid();
+ const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID())->gridView().level();
+
+ HierarchicLevelIteratorType endHierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
+ HierarchicLevelIteratorType hierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
+ for (; hierIt!=endHierIt; ++hierIt) {
+ const Element& fineElem = *hierIt;
+ addLocalDofs(coarseElem, fineElem, coarsePatch, patchDofs);
+ }
+ }
+
+ void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
+ // insert local dofs
+ const auto& lfe = cache_.get(fineElem.type());
+ const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID())->gridView();
+ const auto& indexSet = gridView.indexSet();
+
+ for (size_t i=0; i<lfe.localBasis().size(); ++i) {
+ auto localIdx = lfe.localCoefficients().localKey(i).subEntity();
+ int dofIndex = indexSet.subIndex(fineElem, localIdx, dim);
+ patchDofs[dofIndex][0] = false;
+ }
+
+ const auto coarseLevel = coarseElem.element().level();
+
+ // search for parts of the global and inner boundary
+ for (const auto& is : intersections(gridView, fineElem)) {
+ bool isPatchBoundary = false;
+
+ if (is.neighbor()) {
+ auto outsideFather = is.outside();
+
+ while (outsideFather.level()>coarseLevel) {
+ outsideFather = outsideFather.father();
+ }
+
+ BodyElement outside(coarseElem.bodyID(), outsideFather);
+
+ if (!coarsePatch.count(elementIndex(outside))) {
+ isPatchBoundary = true;
+ }
+ } else {
+ isPatchBoundary = true;
+
+ auto outsideFather = is.outside();
+
+ while (outsideFather.level()>coarseLevel) {
+ outsideFather = outsideFather.father();
+ }
+
+ BodyElement outside(coarseElem.bodyID(), outsideFather);
+
+ size_t faceIdx = faceIndex(coarseElem.bodyID(), is);
+
+ if (neighborElements_.count(faceIdx)) {
+ const auto& neighbors = neighborElements_[faceIdx];
+ for (size_t i=0; i<neighbors.size(); i++) {
+ const auto& neighbor = neighbors[i];
+
+ if (coarsePatch.count(elementIndex(outside))) {
+ isPatchBoundary = false;
+ }
+ }
+
+ }
+ }
+
+ if (isPatchBoundary) {
+ const auto& localCoefficients = cache_.get(fineElem.type()).localCoefficients();
+
+ for (size_t i=0; i<localCoefficients.size(); i++) {
+ auto entity = localCoefficients.localKey(i).subEntity();
+ auto codim = localCoefficients.localKey(i).codim();
+
+ if (containsInsideSubentity(fineElem, is, entity, codim)) {
+ int dofIndex = indexSet.subIndex(fineElem, entity, dim);
+ patchDofs[dofIndex][0] = true;
+ }
+ }
+ }
+ }
+ }
+
+
+ template <class IndexSet, class Intersection>
+ void computeIntersectionDofs(const IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs, bool inside = true) {
+ intersectionDofs.clear();
+
+ Element* elem;
+ size_t intersectionIndex;
+
+ if (inside) {
+ elem = &intersection.inside();
+ intersectionIndex = intersection.indexInInside();
+ } else {
+ elem = &intersection.outside();
+ intersectionIndex = intersection.indexInOutside();
+ }
+
+ const int dimElement = Element::dimension;
+ const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem->type());
+
+ for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
+ size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
+ size_t globalIdx = idxSet.subIndex(*elem, idxInElement, dimElement);
+
+ intersectionDofs.insert(globalIdx);
+ }
+ }
+};
+#endif
diff --git a/src/spatial-solving/preconditioners/test/CMakeLists.txt b/src/spatial-solving/preconditioners/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3a5964adeb30007351d7b9a5c14236191cbe0222
--- /dev/null
+++ b/src/spatial-solving/preconditioners/test/CMakeLists.txt
@@ -0,0 +1,12 @@
+# Put your test in here if it needs access to external grids
+set(FAULTNETWORKS_PRECONDITIONERS_TESTS
+ levelfaultpreconditionertest
+ )
+
+set(TESTS ${FAULTNETWORKS_PRECONDITIONERS_TESTS})
+
+# Setup targets, register tests, and add dune flags
+foreach(_test ${TESTS})
+ #dune_add_test(SOURCES ${_test}.cc)
+endforeach()
+
diff --git a/src/spatial-solving/preconditioners/test/Makefile.am b/src/spatial-solving/preconditioners/test/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..a4c19c2c84cac95071500f9fa9c79aa2677057d1
--- /dev/null
+++ b/src/spatial-solving/preconditioners/test/Makefile.am
@@ -0,0 +1,42 @@
+# which tests where program to build and run are equal
+NORMALTESTS = levelfaultpreconditionertest
+
+# list of tests to run
+TESTS = $(NORMALTESTS)
+
+# programs just to build when "make check" is used
+check_PROGRAMS = $(NORMALTESTS)
+
+# define the programs
+levelfaultpreconditionertest_SOURCES = levelfaultpreconditionertest.cc
+
+
+levelfaultpreconditionertest_CPPFLAGS = $(AM_CPPFLAGS) \
+ $(DUNEMPICPPFLAGS) \
+ $(UG_CPPFLAGS)
+
+# The libraries have to be given in reverse order (most basic libraries
+# last). Also, due to some misunderstanding, a lot of libraries include the
+# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
+# here as well.
+levelfaultpreconditionertest_LDADD = \
+ $(DUNE_LDFLAGS) $(DUNE_LIBS) \
+ $(UG_LDFLAGS) $(UG_LIBS) \
+ $(DUNEMPILIBS) \
+ $(LDADD)
+levelfaultpreconditionertest_LDFLAGS = $(AM_LDFLAGS) \
+ $(DUNEMPILDFLAGS) \
+ $(UG_LDFLAGS) \
+ $(DUNE_LDFLAGS)
+
+# don't follow the full GNU-standard
+# we need automake 1.9
+AUTOMAKE_OPTIONS = foreign 1.9
+
+# pass most important options when "make distcheck" is used
+DISTCHECK_CONFIGURE_FLAGS = --with-dune-common=$(DUNE_COMMON_ROOT) --with-dune-geometry=$(DUNE_GEOMETRY_ROOT) --with-dune-grid=$(DUNE_GRID_ROOT) --with-dune-istl=$(DUNE_ISTL_ROOT) --with-dune-localfunctions=$(DUNE_LOCALFUNCTIONS_ROOT) --with-dune-solvers=$(DUNE_SOLVERS_ROOT) --with-dune-fufem=$(DUNE_FUFEM_ROOT) CXX="$(CXX)" CC="$(CC)"
+
+EXTRA_DIST = CMakeLists.txt
+
+include $(top_srcdir)/am/global-rules
+
diff --git a/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.cc b/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..69b2283a34652d187aa57bf7bb0623e0168cb519
--- /dev/null
+++ b/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.cc
@@ -0,0 +1,274 @@
+#include <config.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
+
+#include <dune/faultnetworks/interfacenetwork.hh>
+#include <dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh>
+#include <dune/faultnetworks/preconditioners/levelfaultpreconditioner.hh>
+#include <dune/faultnetworks/preconditioners/multilevelfaultpreconditioner.hh>
+#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
+#include <dune/faultnetworks/faultp1nodalbasis.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+const int dim = 2;
+
+const int coarseResolution = 2;
+const int fineResolution = 3;
+const int exactResolution = 4;
+
+const int faultCount = 3;
+
+
+typedef typename Dune::UGGrid<dim> GridType;
+typedef typename GridType::LevelGridView GridView;
+typedef typename GridType::LevelIntersection Intersection;
+typedef typename GridType::template Codim<0>::Entity Element;
+
+static const int dimElement = Element::dimension;
+
+void gridLevelOutput(const GridType& grid, const int level) {
+ std::cout << "Grid elements on lvl " << level << std::endl << std::endl;
+
+ typedef typename GridView::template Codim<0>::Iterator ElementLevelIterator;
+ GridView gridView(grid.levelGridView(level));
+
+ //Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > mapper(grid, level);
+ Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGVertexLayout > vMapper(grid, level);
+
+ Dune::BitSetVector<1> vertexPrinted(vMapper.size(), false);
+
+ ElementLevelIterator elemIt = gridView.template begin<0>();
+ ElementLevelIterator elemEndIt = gridView.template end<0>();
+ for (; elemIt!=elemEndIt; ++elemIt) {
+ const Element& elem = *elemIt;
+ const auto& geometry = elem.geometry();
+
+ for(int i=0; i<geometry.corners(); ++i) {
+ size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
+ if (!vertexPrinted[globalIdx][0]) {
+ const auto& vertex = geometry.corner(i);
+ print(vertex, "elem corner:" + std::to_string(globalIdx));
+ vertexPrinted[globalIdx][0] = true;
+ }
+ }
+ }
+ std::cout << "----------------" << std::endl << std::endl;
+
+ /*
+ ElementLevelIterator elemIt = gridView.template begin<0>();
+ ElementLevelIterator elemEndIt = gridView.template end<0>();
+ for (; elemIt!=elemEndIt; ++elemIt) {
+ const Element& elem = *elemIt;
+ const auto& geometry = elem.geometry();
+
+ std::cout << "Element idx: " << mapper.index(elem) << std::endl;
+
+ for(int i=0; i<geometry.corners(); ++i) {
+ size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
+ const auto& vertex = geometry.corner(i);
+ print(vertex, "elem corner:" + std::to_string(globalIdx));
+ }
+ std::cout << "----------------" << std::endl << std::endl;
+ }*/
+
+}
+
+
+void faultInterfaceOutput(const FaultInterface<GridView>& faultInterface) {
+ std::cout << "Level: " << faultInterface.level() << " global faultID: " << faultInterface.id() << std::endl;
+
+ const std::set<size_t>& interfaceDofs = faultInterface.getInterfaceDofs();
+ print(interfaceDofs, "");
+}
+
+void interfaceNetworkLevelOutput(const InterfaceNetwork<GridType>& interfaceNetwork, const int level) {
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
+
+ std::cout << "LevelInterfaceNetwork level: " << level << std::endl;
+ for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
+ std::cout << "local faultID : " << i << " ";
+ faultInterfaceOutput(*(levelInterfaceNetwork.getInterface(i)));
+ }
+
+ std::cout << "------------------------------" << std::endl << std::endl;
+}
+
+void computeIntersectionDofs(const GridView& gridView, const Intersection& intersection, std::set<size_t>& intersectionDofs) {
+ intersectionDofs.clear();
+
+ // loop over all vertices of the intersection
+ const Element& insideElement = intersection.inside();
+
+ const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
+ for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
+ size_t globalIdx = gridView.indexSet().subIndex(insideElement, idxInElement, dimElement);
+ intersectionDofs.insert(globalIdx);
+ }
+}
+
+void prolongFaultDofs(const GridType& grid, const FaultInterface<GridView>& faultInterface, const int toLevel, std::set<size_t>& prolongedDofs) {
+ prolongedDofs.clear();
+
+ const std::vector<Intersection>& faces = faultInterface.faces();
+ const GridView gridView(grid.levelGridView(toLevel));
+
+ for (size_t i=0; i<faces.size(); i++) {
+ const Intersection& intersection = faces[i];
+ Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
+
+ typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(toLevel);
+ typename Face<GridType>::HierarchicIterator exactEnd = face.hend(toLevel);
+ for(; exactIt!=exactEnd; ++exactIt) {
+ if (exactIt->level() == toLevel) {
+ std::set<size_t> faceDofs;
+ computeIntersectionDofs(gridView, *(exactIt->intersection()), faceDofs);
+
+ prolongedDofs.insert(faceDofs.begin(), faceDofs.end());
+ }
+ }
+ }
+}
+
+void prolongInterfaceNetworkLevelDofs(const GridType& grid, const InterfaceNetwork<GridType>& interfaceNetwork, const int level, const int toLevel, std::vector<std::set<size_t>>& prolongedFaultDofs) {
+ prolongedFaultDofs.resize(0);
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
+
+ std::cout << "Correct prolonged LevelInterfaceNetwork level " << level << " to level " << toLevel << std::endl;
+ for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
+ std::set<size_t> prolongedDofs;
+ prolongFaultDofs(grid, *(levelInterfaceNetwork.getInterface(i)), toLevel, prolongedDofs);
+
+ print(prolongedDofs, "local faultID: " + std::to_string(i));
+ prolongedFaultDofs.push_back(prolongedDofs);
+ }
+
+ std::cout << "------------------------------" << std::endl << std::endl;
+}
+
+int main(int argc, char** argv) try
+{
+ typedef double ctype;
+
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::ofstream out("/home/mi/podlesny/software/dune/faultnetworks/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.log");
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt!
+
+ const int fineLevelIdx = fineResolution - coarseResolution;
+ const int exactLevelIdx = exactResolution - coarseResolution;
+
+ std::vector<int> faultToLevel(faultCount);
+ for (size_t i=0; i<faultToLevel.size(); ) {
+ for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
+ faultToLevel[i] = j;
+ i++;
+ }
+
+ for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
+ faultToLevel[i] = j;
+ i++;
+ }
+ }
+
+ print(faultToLevel, "faultToLevel: ");
+
+ MultilevelUniformFaultFactory<GridType> faultFactory(faultCount, faultToLevel, coarseResolution, exactLevelIdx);
+ const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
+
+ /*
+ const GridType& grid = faultFactory.grid();
+
+ gridLevelOutput(grid, 0);
+ gridLevelOutput(grid, fineLevelIdx);
+ gridLevelOutput(grid, exactLevelIdx);
+
+ interfaceNetworkLevelOutput(interfaceNetwork, 0);
+ interfaceNetworkLevelOutput(interfaceNetwork, fineLevelIdx);
+ interfaceNetworkLevelOutput(interfaceNetwork, exactLevelIdx); */
+
+ typedef typename Dune::BlockVector<Dune::FieldVector<ctype,1>> VectorType;
+ typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<ctype, 1, 1>> MatrixType;
+
+ typedef FaultP1NodalBasis<GridView, ctype> DGBasis;
+ typedef typename DGBasis::LocalFiniteElement DGFE;
+
+ typedef IntersectionJumpMassAssembler<GridView, DGFE, DGFE> LocalInterfaceAssembler;
+ typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
+ OscMatrixType B(2,2);
+ B[0][0] = 1;
+ B[1][1] = 1;
+ B[0][1] = -1;
+ B[1][0] = -1;
+
+ //--------------------------------------------------------------
+ // Multilevelfaultpreconditioner
+ //--------------------------------------------------------------
+
+ // setting local assemblers for laplace operator
+ LocalInterfaceAssembler localInterfaceAssembler(B);
+ typedef LaplaceAssembler<GridType, DGFE, DGFE> LocalLaplaceAssembler;
+ LocalLaplaceAssembler localAssembler;
+
+ Dune::BitSetVector<1> activeLevels(fineLevelIdx+1, true);
+ print(activeLevels, "activeLevels: ");
+
+ std::vector<std::shared_ptr<LocalLaplaceAssembler>> localAssemblers;
+ std::vector<std::shared_ptr<LocalInterfaceAssembler>> localIntersectionAssemblers;
+
+ localAssemblers.resize(activeLevels.size());
+ localIntersectionAssemblers.resize(activeLevels.size());
+
+ for (size_t i=0; i<activeLevels.size(); i++) {
+ localAssemblers[i] = std::make_shared<LocalLaplaceAssembler>(localAssembler);
+ localIntersectionAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(localInterfaceAssembler);
+ }
+
+ MultilevelFaultPreconditioner<DGBasis, LocalLaplaceAssembler, LocalInterfaceAssembler, MatrixType, VectorType> preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers);
+
+ MatrixType matrix;
+ VectorType mgX(preconditioner.basis().size());
+ VectorType mgRhs(preconditioner.basis().size());
+ mgRhs = 0;
+
+ preconditioner.setProblem(matrix, mgX, mgRhs);
+ preconditioner.iterate();
+
+ const VectorType& mgSol = preconditioner.getSol();
+ print(mgSol, "mgSol: ");
+
+ step(0);
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(1);
+ LevelFaultPreconditioner<DGBasis, LocalLaplaceAssembler, LocalInterfaceAssembler, MatrixType, VectorType> levelFaultPreconditioner(levelInterfaceNetwork, localAssembler, localInterfaceAssembler);
+
+ step(1);
+
+ VectorType x(levelFaultPreconditioner.basis().size());
+ VectorType rhs(levelFaultPreconditioner.basis().size());
+ rhs = 0;
+
+ step(2);
+
+
+ levelFaultPreconditioner.setProblem(matrix, x, rhs);
+
+ step(3);
+
+ levelFaultPreconditioner.iterate();
+
+ const VectorType& sol = levelFaultPreconditioner.getSol();
+
+ print(sol, "sol: ");
+
+ std::cout.rdbuf(coutbuf);
+} catch (Dune::Exception e){
+ std::cout << e << std::endl;
+}
diff --git a/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.log b/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.log
new file mode 100644
index 0000000000000000000000000000000000000000..a66d14bfdf5e4aed77911c4071d9eda95e757075
--- /dev/null
+++ b/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.log
@@ -0,0 +1,26 @@
+faultToLevel:
+1 0 1
+
+Coarse and fine grid were generated!
+activeLevels:
+1 1
+
+Initializing local fine grid corrections!
+
+Total setup time: 0.000599407 seconds.
+Initializing local fine grid corrections!
+
+Total setup time: 0.00250133 seconds.
+mgSol:
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+
+ Step 0!
+Initializing local fine grid corrections!
+
+Total setup time: 0.00260553 seconds.
+ Step 1!
+ Step 2!
+ Step 3!
+sol:
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index b4a5efd5ffe2b2aa9496da24892412329dea1e62..607ac3658ee0c185cd82f5a949a3dc5e97164af3 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -17,12 +17,13 @@
#include <dune/tnnmg/projections/obstacledefectprojection.hh>
#include <dune/tnnmg/linearsolvers/fastamgsolver.hh>
-//#include "tnnmg/functional.hh"
#include "tnnmg/tnnmgstep.hh"
#include "tnnmg/linearization.hh"
#include "tnnmg/linesearchsolver.hh"
#include "tnnmg/localbisectionsolver.hh"
+//#include "preconditioners/supportpatchfactory.hh"
+
template <class FunctionalTEMPLATE, class BitVectorType>
class SolverFactory {
public:
diff --git a/src/spatial-solving/solverfactory_old.cc b/src/spatial-solving/solverfactory_old.cc
deleted file mode 100644
index cbbe5ef1d7cbabc1c2ecd56677dad245192ed600..0000000000000000000000000000000000000000
--- a/src/spatial-solving/solverfactory_old.cc
+++ /dev/null
@@ -1,78 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/solvers/solvers/solver.hh>
-
-#include "solverfactory_old.hh"
-
-
-template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
-SolverFactoryOld<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactoryOld(
- const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
- const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes)
- : nBodyAssembler_(nBodyAssembler),
- baseEnergyNorm_(baseSolverStep_),
- baseSolver_(&baseSolverStep_,
- parset.get<size_t>("linear.maximumIterations"),
- parset.get<double>("linear.tolerance"), &baseEnergyNorm_,
- Solver::QUIET),
- transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
-
- /* baseSolverStep_.obstacles_ = using HasObstacle = Dune::BitSetVector<BlockSize>;
- using Obstacle = BoxConstraint<Field, BlockSize>;
- const HasObstacle* hasObstacle_;
- const std::vector<Obstacle>* obstacles_;*/
-
- multigridStep_ = std::make_shared<Step>();
- // tnnmg iteration step
- multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
- multigridStep_->setIgnore(ignoreNodes);
- multigridStep_->setBaseSolver(baseSolver_);
- multigridStep_->setSmoother(&presmoother_, &postsmoother_);
- multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
- multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
- multigridStep_->setVerbosity(NumProc::QUIET);
-
- // create the transfer operators
- const int nCouplings = nBodyAssembler_.nCouplings();
- const auto grids = nBodyAssembler_.getGrids();
-
- for (size_t i=0; i<transferOperators_.size(); i++) {
- transferOperators_[i] = std::make_shared<TransferOperator>();
- }
-
- std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
- std::vector<const MatrixType*> mortarTransfer(nCouplings);
- std::vector<std::array<int,2> > gridIdx(nCouplings);
-
- for (int j=0; j<nCouplings; j++) {
- fineHasObstacle[j] = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary().getVertices();
- mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j]->mortarLagrangeMatrix();
- gridIdx[j] = nBodyAssembler_.coupling_[j].gridIdx_;
- }
-
- TransferOperator::setupHierarchy(transferOperators_,
- grids,
- mortarTransfer,
- nBodyAssembler_.localCoordSystems_,
- fineHasObstacle,
- gridIdx);
-
- multigridStep_->setTransferOperators(transferOperators_);
-}
-
-template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
-auto SolverFactoryOld<DeformedGrid, Nonlinearity, MatrixType, VectorType>::getStep()
- -> std::shared_ptr<Step> {
- return multigridStep_;
-}
-
-#include "solverfactory_tmpl_old.cc"
diff --git a/src/spatial-solving/solverfactory_old.hh b/src/spatial-solving/solverfactory_old.hh
deleted file mode 100644
index 1c383419e1ac32eced29512d04260881f7157a82..0000000000000000000000000000000000000000
--- a/src/spatial-solving/solverfactory_old.hh
+++ /dev/null
@@ -1,66 +0,0 @@
-#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORYOLD_HH
-#define SRC_SPATIAL_SOLVING_SOLVERFACTORYOLD_HH
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
-#include <dune/solvers/iterationsteps/blockgsstep.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-
-//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
-
-#include <dune/contact/assemblers/nbodyassembler.hh>
-
-//#include <dune/contact/estimators/hierarchiccontactestimator.hh>
-#include <dune/contact/solvers/nsnewtonmgstep.hh>
-#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
-#include <dune/contact/solvers/contactobsrestrict.hh>
-#include <dune/contact/solvers/contacttransferoperatorassembler.hh>
-#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
-
-#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
-
-template <class DeformedGridTEMPLATE, class NonlinearityTEMPLATE, class MatrixType, class VectorType>
-class SolverFactoryOld {
-//protected:
-// const size_t dim = DeformedGrid::dimension;
-
-public:
- using Matrix = MatrixType;
- using Vector = VectorType;
- using DeformedGrid = DeformedGridTEMPLATE;
- using Nonlinearity = NonlinearityTEMPLATE;
-
-public:
- using Step = Dune::Contact::NonSmoothNewtonMGStep<MatrixType, VectorType>;
- using TransferOperator = Dune::Contact::NonSmoothNewtonContactTransfer<VectorType>;
-
- SolverFactoryOld(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
- const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes);
-
- std::shared_ptr<Step> getStep();
-
- const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& getNBodyAssembler() const {
- return nBodyAssembler_;
- }
-
-private:
- const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler_;
-
- //ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
- BlockGSStep<MatrixType, VectorType> baseSolverStep_;
- EnergyNorm<MatrixType, VectorType> baseEnergyNorm_;
- LoopSolver<VectorType> baseSolver_;
-
- ProjectedBlockGSStep<MatrixType, VectorType> presmoother_;
- ProjectedBlockGSStep<MatrixType, VectorType> postsmoother_;
- std::shared_ptr<Step> multigridStep_;
-
- std::vector<std::shared_ptr<TransferOperator>> transferOperators_;
-};
-#endif
diff --git a/src/spatial-solving/solverfactory_tmpl_old.cc b/src/spatial-solving/solverfactory_tmpl_old.cc
deleted file mode 100644
index 600132d31b71a6fda13c7d955374bf98d01f9412..0000000000000000000000000000000000000000
--- a/src/spatial-solving/solverfactory_tmpl_old.cc
+++ /dev/null
@@ -1,16 +0,0 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-
-#include <dune/tectonic/globalfriction.hh>
-
-template class SolverFactoryOld<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
-
-/*template class SolverFactory<
- MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;*/
diff --git a/src/spatial-solving/tnnmg/CMakeLists.txt b/src/spatial-solving/tnnmg/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..25682fc9e6c4962dd1d531ad1e5917944aca91d3
--- /dev/null
+++ b/src/spatial-solving/tnnmg/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_custom_target(dune-tectonic_spatial-solving_tnnmg_src SOURCES
+ functional.hh
+ linearization.hh
+ linesearchsolver.hh
+ localbisectionsolver.hh
+ localproblem.hh
+ tnnmgstep.hh
+ zerononlinearity.hh
+)
diff --git a/src/spatial-solving/tnnmg/levelpatchpreconditioner.hh b/src/spatial-solving/tnnmg/levelpatchpreconditioner.hh
deleted file mode 100644
index bbe464fdec9c3d9db5c218cecdbe4c448cd06fad..0000000000000000000000000000000000000000
--- a/src/spatial-solving/tnnmg/levelpatchpreconditioner.hh
+++ /dev/null
@@ -1,261 +0,0 @@
-#ifndef LEVEL_PATCH_PRECONDITIONER_HH
-#define LEVEL_PATCH_PRECONDITIONER_HH
-
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-
-#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
-#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
-#include <dune/faultnetworks/localproblem.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-
-template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
-class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- enum Mode {additive, multiplicative};
- enum BoundaryMode {homogeneous, fromIterate};
-
-private:
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
- const BasisType& patchLevelBasis_;
- const LocalAssembler& localAssembler_;
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
- const Mode mode_;
-
- const GridType& grid_;
- const int level_;
- const BasisType basis_;
-
- size_t patchDepth_;
- BoundaryMode boundaryMode_;
-
- MatrixType matrix_;
- std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
-
-public:
-
- // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
- LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
- const BasisType& patchLevelBasis,
- const LocalAssembler& localAssembler,
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
- const Mode mode = LevelPatchPreconditioner::Mode::additive) :
- levelInterfaceNetwork_(levelInterfaceNetwork),
- patchLevelBasis_(patchLevelBasis),
- localAssembler_(localAssembler),
- localInterfaceAssemblers_(localInterfaceAssemblers),
- mode_(mode),
- grid_(levelInterfaceNetwork_.grid()),
- level_(levelInterfaceNetwork_.level()),
- basis_(levelInterfaceNetwork_)
- {
- setPatchDepth();
- setBoundaryMode();
-
- assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
- }
-
- void build() {
- // assemble stiffness matrix for entire level including all faults
- GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
- globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
-
- // set boundary conditions
- Dune::BitSetVector<1> globalBoundaryDofs;
- BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
- constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
- if(!globalBoundaryDofs[i][0])
- continue;
-
- RowType& row = matrix_[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- row[cIt.index()] = 0;
- }
- row[i] = 1;
- }
-
- // init vertexInElements
- const int dim = GridType::dimension;
- typedef typename GridType::template Codim<0>::Entity EntityType;
- std::vector<std::vector<EntityType>> vertexInElements;
-
- const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
- vertexInElements.resize(patchLevelGridView.size(dim));
- for (size_t i=0; i<vertexInElements.size(); i++) {
- vertexInElements[i].resize(0);
- }
-
- typedef typename BasisType::LocalFiniteElement FE;
- typedef typename GridView::template Codim <0>::Iterator ElementLevelIterator;
- ElementLevelIterator endElemIt = patchLevelGridView.template end <0>();
- for (ElementLevelIterator it = patchLevelGridView. template begin <0>(); it!=endElemIt; ++it) {
-
- // compute coarseGrid vertexInElements
- const FE& coarseFE = patchLevelBasis_.getLocalFiniteElement(*it);
-
- for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
- int dofIndex = patchLevelBasis_.indexInGridView(*it, i);
- vertexInElements[dofIndex].push_back(*it);
- }
- }
-
- localProblems_.resize(vertexInElements.size());
-
- std::cout << std::endl;
- std::cout << "---------------------------------------------" << std::endl;
- std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
-
- // init local fine level corrections
- Dune::Timer timer;
- timer.reset();
- timer.start();
-
- const int patchLevel = patchLevelBasis_.faultNetwork().level();
- for (size_t i=0; i<vertexInElements.size(); i++) {
- //std::cout << i << std::endl;
- //std::cout << "---------------" << std::endl;
-
- SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, basis_, patchLevel, level_, vertexInElements, i, patchDepth_);
- std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
- VectorType rhs;
- rhs.resize(matrix_.N());
- rhs = 0;
-
- //print(localToGlobal, "localToGlobal: ");
- //print(boundaryDofs, "boundaryDofs: ");
-
- localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);
-
- /*
- if ((i+1) % 10 == 0) {
- std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
- }*/
- }
-
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- std::cout << "---------------------------------------------" << std::endl;
- timer.stop();
- }
-
- void setPatchDepth(const size_t patchDepth = 0) {
- patchDepth_ = patchDepth;
- }
-
- void setBoundaryMode(const BoundaryMode boundaryMode = LevelPatchPreconditioner::BoundaryMode::homogeneous) {
- boundaryMode_ = boundaryMode;
- }
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- if (boundaryMode_ == BoundaryMode::homogeneous)
- localProblems_[i]->updateRhs(rhs);
- else
- localProblems_[i]->updateRhsAndBoundary(rhs, x);
- }
- }
-
- virtual void iterate() {
- if (mode_ == additive)
- iterateAdd();
- else
- iterateMult();
- }
-
- void iterateAdd() {
- *(this->x_) = 0;
-
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->solve(it);
- localProblems_[i]->prolong(it, x);
-
- /*if (i==5) {
- writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
- }*/
-
- *(this->x_) += x;
- }
- }
-
- void iterateMult() {
- *(this->x_) = 0;
-
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
- VectorType updatedRhs(*(this->rhs_));
- matrix_.mmv(*(this->x_), updatedRhs);
-
- //step(i);
- //print(updatedRhs, "localRhs: ");
- //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));
-
- if (boundaryMode_ == BoundaryMode::homogeneous)
- localProblems_[i]->updateRhs(updatedRhs);
- else
- localProblems_[i]->updateRhsAndBoundary(updatedRhs, *(this->x_));
-
- localProblems_[i]->solve(it);
- localProblems_[i]->prolong(it, x);
-
-
- //print(it, "it: ");
- //print(x, "x: ");
-
- //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));
-
- /*if (i==5) {
- writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
- }*/
-
- *(this->x_) += x;
- }
- }
-
- const BasisType& basis() const {
- return basis_;
- }
-
- const GridView& gridView() const {
- return basis_.getGridView();
- }
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
- return levelInterfaceNetwork_;
- }
-
- size_t size() const {
- return localProblems_.size();
- }
-};
-
-#endif
-
diff --git a/src/spatial-solving/tnnmg/supportpatchfactory.hh b/src/spatial-solving/tnnmg/supportpatchfactory.hh
deleted file mode 100644
index 20754c0d5c371b5d9a715778a2b9d7e87cace00d..0000000000000000000000000000000000000000
--- a/src/spatial-solving/tnnmg/supportpatchfactory.hh
+++ /dev/null
@@ -1,253 +0,0 @@
-#ifndef SUPPORT_PATCH_FACTORY_HH
-#define SUPPORT_PATCH_FACTORY_HH
-
-#include<queue>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/fufem/referenceelementhelper.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-
-template <class BasisType>
-class SupportPatchFactory
-{
- protected:
- typedef typename BasisType::GridView GV;
- typedef typename BasisType::LocalFiniteElement LFE;
-
- typedef typename GV::Grid GridType;
- typedef typename GridType::ctype ctype;
- static const int dim = GridType::dimension;
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::template Codim<dim>::Entity Vertex;
-
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper;
-
- private:
- const BasisType& coarseBasis_;
- const BasisType& fineBasis_;
-
- const int coarseLevel_;
- const int fineLevel_;
-
- const std::vector< std::vector <Element>>& vertexInElements_;
- const int centerVertexIdx_;
- const int patchDepth_;
-
-
- const GridType& grid;
- const GV& coarseGridView;
- const GV& fineGridView;
-
- std::vector<int> localToGlobal;
- Dune::BitSetVector<1> boundaryDofs;
- std::vector<Element> fineRegionElements;
-
- ElementMapper coarseMapper;
-
- public:
- //setup
- SupportPatchFactory(const BasisType& coarseBasis, const BasisType& fineBasis, const int coarseLevel, const int fineLevel, const std::vector< std::vector <Element>>& vertexInElements, const int centerVertexIdx, const int patchDepth) :
- coarseBasis_(coarseBasis),
- fineBasis_(fineBasis),
- coarseLevel_(coarseLevel),
- fineLevel_(fineLevel),
- vertexInElements_(vertexInElements),
- centerVertexIdx_(centerVertexIdx),
- patchDepth_(patchDepth),
- grid(coarseBasis_.getGridView().grid()),
- coarseGridView(coarseBasis_.getGridView()),
- fineGridView(fineBasis_.getGridView()),
- coarseMapper(grid, coarseLevel_)
- {
- fineRegionElements.resize(0);
-
- std::vector<Element> coarseElements;
- Dune::BitSetVector<1> visited(coarseMapper.size());
- Dune::BitSetVector<1> vertexVisited(vertexInElements_.size());
- visited.unsetAll();
- vertexVisited.unsetAll();
-
- std::set<int> localDofs;
- std::set<int> localBoundaryDofs;
-
- addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
-
- // construct coarse patch
- for (int depth=1; depth<patchDepth_; depth++) {
- std::vector<Element> coarseElementNeighbors;
-
- for (size_t i=0; i<coarseElements.size(); ++i) {
- const Element& elem = coarseElements[i];
- const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
-
- for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
- int dofIndex = coarseBasis_.indexInGridView(elem, j);
-
- if (!vertexVisited[dofIndex][0]) {
- addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
- }
- }
- }
-
- coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
- }
-
-
- // construct fine patch
- for (size_t i=0; i<coarseElements.size(); ++i) {
- const Element& elem = coarseElements[i];
- addFinePatchElements(elem, visited, localDofs, localBoundaryDofs);
- }
-
-
-
- localToGlobal.resize(localDofs.size());
- boundaryDofs.resize(localDofs.size());
- boundaryDofs.unsetAll();
-
- std::set<int>::iterator dofIt = localDofs.begin();
- std::set<int>::iterator dofEndIt = localDofs.end();
- size_t i=0;
- for (; dofIt != dofEndIt; ++dofIt) {
- int localDof = *dofIt;
- localToGlobal[i] = localDof;
-
- if (localBoundaryDofs.count(localDof)) {
- boundaryDofs[i][0] = true;
- }
- i++;
- }
- }
-
- std::vector<int>& getLocalToGlobal() {
- return localToGlobal;
- }
-
- std::vector<Element>& getRegionElements() {
- return fineRegionElements;
- }
-
- Dune::BitSetVector<1>& getBoundaryDofs() {
- return boundaryDofs;
- }
-
-
-private:
-
-
- void print(const Dune::BitSetVector<1>& x, std::string message){
- std::cout << message << std::endl;
- for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i][0] << " ";
- }
- std::cout << std::endl << std::endl;
- }
-
- void addVertexSupport(const int vertexIdx, std::vector<Element>& coarseElements, Dune::BitSetVector<1>& visited, Dune::BitSetVector<1>& vertexVisited) {
- const std::vector<Element>& vertexInElement = vertexInElements_[vertexIdx];
-
- for (size_t i=0; i<vertexInElement.size(); i++) {
- const Element& elem = vertexInElement[i];
- const int elemIdx = coarseMapper.index(elem);
-
- if (!visited[elemIdx][0]) {
- coarseElements.push_back(elem);
- visited[elemIdx][0] = true;
- }
- }
- vertexVisited[vertexIdx][0] = true;
- }
-
- bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
- return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
- }
-
-
- void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
- HierarchicLevelIteratorType endHierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
- for (HierarchicLevelIteratorType hierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
-
- const Element& fineElem = *hierIt;
- fineRegionElements.push_back(fineElem);
-
- addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs);
- }
- }
-
- void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
- const LFE& fineLFE = fineBasis_.getLocalFiniteElement(fineElem);
-
- /*
- const auto& fineGeometry = fineElem.geometry();
- const auto& coarseGeometry = coarseElem.geometry();
-
- const auto& ref = Dune::ReferenceElements<ctype, dim>::general(fineElem.type()); */
-
- // insert local dofs
- for (size_t i=0; i<fineLFE.localBasis().size(); ++i) {
- int dofIndex = fineBasis_.index(fineElem, i);
- localDofs.insert(dofIndex);
- }
-
-
- // search for parts of the global and inner boundary
- typename GridType::LevelIntersectionIterator neighborIt = fineGridView.ibegin(fineElem);
- typename GridType::LevelIntersectionIterator neighborItEnd = fineGridView.iend(fineElem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const typename GridType::LevelIntersection& intersection = *neighborIt;
-
- bool isInnerBoundary = false;
- if (intersection.neighbor()) {
- Element outsideFather = intersection.outside();
-
- while (outsideFather.level()>coarseLevel_) {
- outsideFather = outsideFather.father();
- }
-
- if (!inCoarsePatch[coarseMapper.index(outsideFather)][0]) {
- isInnerBoundary = true;
- }
- }
-
- if (intersection.boundary() or isInnerBoundary) {
- typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
- const LocalCoefficients* localCoefficients = &fineBasis_.getLocalFiniteElement(intersection.inside()).localCoefficients();
-
- for (size_t i=0; i<localCoefficients->size(); i++) {
- unsigned int entity = localCoefficients->localKey(i).subEntity();
- unsigned int codim = localCoefficients->localKey(i).codim();
-
- if (containsInsideSubentity(fineElem, intersection, entity, codim))
- localBoundaryDofs.insert(fineBasis_.index(intersection.inside(), i));
- }
- }
- }
-
- /* add interior coarse level vertices to boundary dofs
- for(int i=0; i<coarseGeometry.corners(); ++i) {
- const auto& local = fineGeometry.local(coarseGeometry.corner(i));
-
- if (!ref.checkInside(local))
- continue;
-
- const int localBasisSize = fineLFE.localBasis().size();
- std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
- fineLFE.localBasis().evaluateFunction(local, localBasisRep);
-
- for(int k=0; k<localBasisSize; ++k) {
- if (localBasisRep[k]==1) {
- int dofIndex = fineBasis_.index(fineElem, k);
- localBoundaryDofs.insert(dofIndex);
- break;
- }
- }
- } */
-
- }
-};
-#endif
diff --git a/src/tests/gridgluefrictiontest.cc b/src/tests/gridgluefrictiontest.cc
index c33a84937204f1649228608d9fb08f640e034fc9..d31f58d3490ba44c1012a520af030400f1787ef3 100644
--- a/src/tests/gridgluefrictiontest.cc
+++ b/src/tests/gridgluefrictiontest.cc
@@ -399,6 +399,43 @@ int main(int argc, char *argv[]) { try {
std::cout << "---------"<< std::endl;
}*/
+ for (const auto& elem : elements(gridView0)) {
+ std::cout << "seed element corners:" << std::endl;
+
+ const auto& eRefElement = Dune::ReferenceElements<double, dim>::general(elem.type());
+ size_t vertices = eRefElement.size(dim);
+
+ for (size_t j=0; j<vertices; j++) {
+ print(elem.geometry().corner(j), "corner: ");
+ }
+
+
+ for (auto&& is : intersections(gridView0, elem)) {
+ const auto& isGeo = is.geometry();
+
+ if (!is.neighbor())
+ continue;
+
+ std::cout << "neighbor corners:" << std::endl;
+
+ const auto& outside = is.outside();
+
+ const auto& refElement = Dune::ReferenceElements<double, dim>::general(outside.type());
+ size_t nVertices = refElement.size(dim);
+
+ const auto& isRefElement = Dune::ReferenceElements<double, dim-1>::general(isGeo.type());
+
+ for (size_t j=0; j<nVertices; j++) {
+ print(outside.geometry().corner(j), "corner: ");
+ const auto& local = elem.geometry().local(outside.geometry().corner(j));
+
+ print(local, "local:");
+ }
+
+ }
+
+ break;
+ }
bool passed = true;