diff --git a/src/obstaclesolver.hh b/src/obstaclesolver.hh
deleted file mode 100644
index 49b4d7cb2c4c1d75055aec2bc12a22d28a1535e0..0000000000000000000000000000000000000000
--- a/src/obstaclesolver.hh
+++ /dev/null
@@ -1,48 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_TECTONIC_OBSTACLESOLVER_HH
-#define DUNE_TECTONIC_OBSTACLESOLVER_HH
-
-#include "spatial-solving/tnnmg/localproblem.hh"
-
-
-/**
- * \brief A local solver for quadratic obstacle problems
- *
- * \todo Add concept check for the function interface
- */
-class ObstacleSolver
-{
-public:
- template<class Vector, class Functional, class BitVector>
- constexpr void operator()(Vector& x, const Functional& f, const BitVector& ignore) const
- {
- auto& A = f.quadraticPart();
- auto& r = f.linearPart();
-
- LocalProblem<std::decay_t<decltype(A)> , std::decay_t<decltype(r)>> localProblem(A, r, ignore);
- Vector newR;
- localProblem.getLocalRhs(x, newR);
-
- auto directSolver = Dune::Solvers::UMFPackSolver<std::decay_t<decltype(A)>, std::decay_t<decltype(r)>>();
-
- directSolver->setProblem(localProblem.getMat(), x, newR);
- directSolver->preprocess();
- directSolver->solve();
-
- const auto& lower = f.lowerObstacle();
- const auto& upper = f.upperObstacle();
-
- for (size_t i=0; i<x.size(); i++) {
- if (ignore[i][0])
- continue;
-
- if (x[i] < lower[i])
- x[i] = lower[i];
- if (x[i] > upper[i])
- x[i] = upper[i];
- }
- }
-};
-
-#endif
diff --git a/src/solverfactorytest.cc b/src/solverfactorytest.cc
index bbd7aeb09015b3d768a5724d391786d0b630805b..ef3ad5a509aabfe8c60df072fc3b5185cb3529f7 100644
--- a/src/solverfactorytest.cc
+++ b/src/solverfactorytest.cc
@@ -40,10 +40,8 @@
#include "io/vtk.hh"
-#include "spatial-solving/tnnmg/tnnmgstep.hh"
#include "spatial-solving/tnnmg/linesearchsolver.hh"
#include "spatial-solving/preconditioners/nbodycontacttransfer.hh"
-#include "obstaclesolver.hh"
#include "factories/stackedblocksfactory.hh"
diff --git a/src/spatial-solving/CMakeLists.txt b/src/spatial-solving/CMakeLists.txt
index 1970afe64c1da2c2152ba9de21234515d17472fc..e683042b12f4a2e84ed1ea781d911b427d878e63 100644
--- a/src/spatial-solving/CMakeLists.txt
+++ b/src/spatial-solving/CMakeLists.txt
@@ -1,3 +1,2 @@
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
deleted file mode 100644
index 05ec76c58ad389b8aac6c0368f6f551825f2ccb3..0000000000000000000000000000000000000000
--- a/src/spatial-solving/contact/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-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
deleted file mode 100644
index 96ca666c521a49612d334915f26adf359af422fa..0000000000000000000000000000000000000000
--- a/src/spatial-solving/contact/dualmortarcoupling.cc
+++ /dev/null
@@ -1,333 +0,0 @@
-// -*- 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
deleted file mode 100644
index 81c6e628ecd8146abe3847818843aff9a78d9500..0000000000000000000000000000000000000000
--- a/src/spatial-solving/contact/dualmortarcoupling.hh
+++ /dev/null
@@ -1,199 +0,0 @@
-// -*- 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
deleted file mode 100644
index d8ee25a45d13da8fdaa6f5810937c09631388b64..0000000000000000000000000000000000000000
--- a/src/spatial-solving/contact/nbodyassembler.cc
+++ /dev/null
@@ -1,519 +0,0 @@
-// -*- 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
deleted file mode 100644
index ae1db365287d1355a56753b509f3d6bc4e717395..0000000000000000000000000000000000000000
--- a/src/spatial-solving/contact/nbodyassembler.hh
+++ /dev/null
@@ -1,200 +0,0 @@
-// -*- 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
index c5f8a1b52872ab6b9de00bebca17b7077454e541..00e132a99faa7bc54408ea5df6ead98b1a9c4cd5 100644
--- a/src/spatial-solving/preconditioners/CMakeLists.txt
+++ b/src/spatial-solving/preconditioners/CMakeLists.txt
@@ -1,7 +1,7 @@
add_custom_target(dune-tectonic_spatial-solving_preconditioners_src SOURCES
hierarchicleveliterator.hh
- levelglobalpreconditioner.hh
levelpatchpreconditioner.hh
+ localproblem.hh
multilevelpatchpreconditioner.hh
nbodycontacttransfer.hh
supportpatchfactory.hh
diff --git a/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh b/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
deleted file mode 100644
index e593c60daee0051737b80ee1e7d3ee558ba4bcad..0000000000000000000000000000000000000000
--- a/src/spatial-solving/preconditioners/levelglobalpreconditioner.hh
+++ /dev/null
@@ -1,93 +0,0 @@
-#ifndef LEVEL_GLOBAL_PRECONDITIONER_HH
-#define LEVEL_GLOBAL_PRECONDITIONER_HH
-
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-
-//#include <dune/istl/superlu.hh>
-#include <dune/istl/umfpack.hh>
-
-#include "../../data-structures/levelcontactnetwork.hh"
-
-
-template <class LevelContactNetwork, class Solver, class MatrixType, class VectorType>
-class LevelGlobalPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-private:
- using Base = LinearIterationStep<MatrixType, VectorType>;
-
- using ctype = typename LevelContactNetwork::ctype;
-
- const LevelContactNetwork& levelContactNetwork_;
- const int level_;
-
- std::shared_ptr<Solver> solver_;
-
-public:
-
- LevelGlobalPreconditioner(const LevelContactNetwork& levelContactNetwork, const int level = 0) :
- levelContactNetwork_(levelContactNetwork),
- level_(level) {
-
- this->verbosity_ = QUIET;
- }
-
- void setSolver(Solver&& patchSolver) {
- solver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
- }
-
- void iterate() override {
- Dune::Timer timer;
-
- if (this->verbosity_ == FULL) {
- timer.reset();
- timer.start();
- }
-
- *(this->x_) = 0;
-
- auto& step = solver_->getIterationStep();
- step.setProblem(this->mat_, this->x_, this->rhs_);
- step.setIgnore(this->ignoreNodes_);
-
- solver_->check();
- solver_->preprocess();
- solver_->solve();
-
- *(this->x_) = solver_->getSol();
-
- if (this->verbosity_ == FULL) {
- std::cout << "LevelGlobalPreconditioner::iterate() Solving global problem took: " << timer.stop() << " seconds" << std::endl;
- }
- }
-
- /*void iterate() override {
- Dune::Timer timer;
- timer.reset();
- timer.start();
-
- this->x_->resize(matrix_.M());
- *(this->x_) = 0;
-
- #if HAVE_SUPERLU
- Dune::InverseOperatorResult res;
- VectorType rhsCopy(rhs_);
-
- Dune::UMFPack<MatrixType> solver(matrix_);
- solver.apply(*(this->x_), rhsCopy, res);
-
- #else
- #error No SuperLU!
- #endif
-
- //std::cout << "LevelGlobalPreconditioner::iterate() Solving global problem took: " << timer.elapsed() << " seconds" << std::endl;
-
- }*/
-
-};
-#endif
-
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index 6135de9c97dab8a8e4d575718cc54f67f1af10ed..a515f6e5b6405f767168c5734995a0ccca108095 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -13,7 +13,6 @@
#include "../../data-structures/levelcontactnetwork.hh"
#include "supportpatchfactory.hh"
-#include "../tnnmg/localproblem.hh"
#include <dune/localfunctions/lagrange/pqkfactory.hh>
@@ -120,27 +119,6 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
if (!vertexVisited[globalIdx][0]) {
vertexVisited[globalIdx][0] = true;
patchFactory_.build(bodyIdx, e, i, patches_[globalIdx], patchDepth_);
-
- /*print(patches_[globalIdx], "patch:");
-
- size_t c = 0;
- for (size_t j=0; j<levelContactNetwork_.nBodies(); j++) {
- const auto& gv = fineContactNetwork_.body(j)->gridView();
-
- printDofLocation(gv);
-
- Dune::BlockVector<Dune::FieldVector<ctype, 1>> patchVec(gv.size(dim));
- for (size_t l=0; l<patchVec.size(); l++) {
- if (patches_[globalIdx][c++][0]) {
- patchVec[l][0] = 1;
- }
- }
-
- //print(patchVec, "patchVec");
-
- // output patch
- writeToVTK(gv, patchVec, "", "level_" + std::to_string(level_) + "_patch_" + std::to_string(globalIdx) + "_body_" + std::to_string(j));
- }*/
}
}
}
@@ -201,27 +179,14 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
patchSolver_->solve();
*(this->x_) += x;
-
-
-
- /*LocalProblem<MatrixType, VectorType> localProblem(*this->mat_, *this->rhs_, ignore);
- Vector newR;
- localProblem.getLocalRhs(x, newR); */
-
- /*print(ignore, "ignore:");
- print(*this->mat_, "Mat:");
- print(localProblem.getMat(), "localMat:");*/
-
- /*patchSolver_->setProblem(localProblem.getMat(), x, newR);
- patchSolver_->preprocess();
- patchSolver_->solve();
-
- *(this->x_) += x;*/
}
}
void iterateMult() {
*(this->x_) = 0;
+
+ DUNE_THROW(Dune::Exception, "Not implemented!");
+
/*
VectorType x = *(this->x_);
for (size_t i=0; i<patches_.size(); i++) {
diff --git a/src/spatial-solving/tnnmg/localproblem.hh b/src/spatial-solving/preconditioners/localproblem.hh
similarity index 100%
rename from src/spatial-solving/tnnmg/localproblem.hh
rename to src/spatial-solving/preconditioners/localproblem.hh
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 9d6d64e6eaf8391c6fe6e88290c3a76140130ee1..b187da5372f084a446c27b99205c71db66d73d4f 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -19,8 +19,7 @@
#include "nbodycontacttransfer.hh"
#include "levelpatchpreconditioner.hh"
-
-#include "../../utils/debugutils.hh"
+#include "localproblem.hh"
template <class ContactNetwork, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
diff --git a/src/spatial-solving/preconditioners/test/CMakeLists.txt b/src/spatial-solving/preconditioners/test/CMakeLists.txt
deleted file mode 100644
index 3a5964adeb30007351d7b9a5c14236191cbe0222..0000000000000000000000000000000000000000
--- a/src/spatial-solving/preconditioners/test/CMakeLists.txt
+++ /dev/null
@@ -1,12 +0,0 @@
-# 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
deleted file mode 100644
index a4c19c2c84cac95071500f9fa9c79aa2677057d1..0000000000000000000000000000000000000000
--- a/src/spatial-solving/preconditioners/test/Makefile.am
+++ /dev/null
@@ -1,42 +0,0 @@
-# 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
deleted file mode 100644
index 69b2283a34652d187aa57bf7bb0623e0168cb519..0000000000000000000000000000000000000000
--- a/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.cc
+++ /dev/null
@@ -1,274 +0,0 @@
-#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
deleted file mode 100644
index a66d14bfdf5e4aed77911c4071d9eda95e757075..0000000000000000000000000000000000000000
--- a/src/spatial-solving/preconditioners/test/levelfaultpreconditionertest.log
+++ /dev/null
@@ -1,26 +0,0 @@
-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/tnnmg/CMakeLists.txt b/src/spatial-solving/tnnmg/CMakeLists.txt
index 698f4e86b7c4c100104d4aa3372ac418b0f7f20e..d32080a39f211c3f91f60516d56274e94808dd29 100644
--- a/src/spatial-solving/tnnmg/CMakeLists.txt
+++ b/src/spatial-solving/tnnmg/CMakeLists.txt
@@ -4,7 +4,5 @@ add_custom_target(dune-tectonic_spatial-solving_tnnmg_src SOURCES
linearcorrection.hh
linesearchsolver.hh
localbisectionsolver.hh
- localproblem.hh
- tnnmgstep.hh
zerononlinearity.hh
)
diff --git a/src/spatial-solving/tnnmg/linearcorrection.hh b/src/spatial-solving/tnnmg/linearcorrection.hh
index 188d580fa41d56f4d32220591fc502d91a1e2e73..5cf2eae90d6404fa6fb3cf821e1f42797fe742dd 100644
--- a/src/spatial-solving/tnnmg/linearcorrection.hh
+++ b/src/spatial-solving/tnnmg/linearcorrection.hh
@@ -11,7 +11,6 @@
#include <dune/solvers/common/resize.hh>
#include <dune/solvers/solvers/umfpacksolver.hh>
-#include "localproblem.hh"
/**
* \brief linear correction step for use by \ref TNNMGStep
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
deleted file mode 100644
index 1d43966459194bd613a121eb2d020869551ad39d..0000000000000000000000000000000000000000
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ /dev/null
@@ -1,286 +0,0 @@
-#ifndef DUNE_TECTONIC_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
-#define DUNE_TECTONIC_TRUNCATED_NONSMOOTH_NEWTON_MULTIGRID_HH
-
-#include <string>
-#include <sstream>
-#include <vector>
-#include <iomanip>
-
-#include <dune/common/timer.hh>
-
-#include <dune/solvers/common/resize.hh>
-#include "dune/solvers/iterationsteps/iterationstep.hh"
-#include "dune/solvers/iterationsteps/lineariterationstep.hh"
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/solvers/linearsolver.hh>
-#include <dune/solvers/common/canignore.hh>
-
-#include <dune/solvers/solvers/umfpacksolver.hh>
-
-#include "localproblem.hh"
-
-#include "linearcorrection.hh"
-//#include <dune/tnnmg/iterationsteps/linearcorrection.hh>
-
-#include <dune/tectonic/../../src/utils/debugutils.hh>
-
-
-/**
- * \brief One iteration of the TNNMG method
- *
- * \tparam F Functional to minimize
- * \tparam BV Bit-vector type for marking ignored components
- */
-template<class F, class BV, class Linearization,
- class DefectProjection,
- class LineSearchSolver>
-class TNNMGStep :
- public IterationStep<typename F::Vector, BV>
-{
- using Base = IterationStep<typename F::Vector, BV>;
-
-public:
-
- using Vector = typename F::Vector;
- using ConstrainedVector = typename Linearization::ConstrainedVector;
- using ConstrainedMatrix = typename Linearization::ConstrainedMatrix;
- using BitVector = typename Base::BitVector;
- using ConstrainedBitVector = typename Linearization::ConstrainedBitVector;
- using Functional = F;
- using IterativeSolver = Dune::Solvers::IterativeSolver< ConstrainedVector, Dune::Solvers::DefaultBitVector_t<ConstrainedVector> >;
- using LinearSolver = Dune::Solvers::LinearSolver< ConstrainedMatrix, ConstrainedVector >;
-
- /** \brief Constructor with an iterative solver object for the linear correction
- * \param iterativeSolver This is a callback used to solve the constrained linearized system
- * \param projection This is a callback used to compute a projection into a defect-admissible set
- * \param lineSolver This is a callback used to minimize a directional restriction of the functional
- * for computing a damping parameter
- */
- TNNMGStep(const Functional& f,
- Vector& x,
- std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother,
- std::shared_ptr<IterativeSolver> iterativeSolver,
- const DefectProjection& projection,
- const LineSearchSolver& lineSolver)
- : Base(x),
- f_(&f),
- nonlinearSmoother_(nonlinearSmoother),
- linearCorrection_(buildLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
- projection_(projection),
- lineSolver_(lineSolver)
- {}
-
- /** \brief Constructor with a linear solver object for the linear correction
- * \param linearSolver This is a callback used to solve the constrained linearized system
- * \param projection This is a callback used to compute a projection into a defect-admissible set
- * \param lineSolver This is a callback used to minimize a directional restriction of the functional
- * for computing a damping parameter
- */
- TNNMGStep(const Functional& f,
- Vector& x,
- std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother,
- std::shared_ptr<LinearSolver> linearSolver,
- const DefectProjection& projection,
- const LineSearchSolver& lineSolver)
- : Base(x),
- f_(&f),
- nonlinearSmoother_(nonlinearSmoother),
- linearCorrection_(buildLinearCorrection(linearSolver)),
- projection_(projection),
- lineSolver_(lineSolver)
- {}
-
- /** \brief Constructor with a LinearIterationStep object for the linear correction
- * \param linearIterationStep This is a callback used to solve the constrained linearized system
- * \param projection This is a callback used to compute a projection into a defect-admissible set
- * \param lineSolver This is a callback used to minimize a directional restriction of the functional
- * for computing a damping parameter
- */
- TNNMGStep(const Functional& f,
- Vector& x,
- std::shared_ptr<Dune::Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
- std::shared_ptr<Dune::Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
- unsigned int noOfLinearIterationSteps,
- const DefectProjection& projection,
- const LineSearchSolver& lineSolver)
- : Base(x),
- f_(&f),
- nonlinearSmoother_(nonlinearSmoother),
- linearCorrection_(buildLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
- projection_(projection),
- lineSolver_(lineSolver)
- {}
-
- using Base::getIterate;
-
- void preprocess() override
- {
- nonlinearSmoother_->setIgnore(this->ignore());
- nonlinearSmoother_->preprocess();
- }
-
- void setPreSmoothingSteps(std::size_t i)
- {
- preSmoothingSteps_ = i;
- }
-
- auto energy(const Vector& v) const {
- Vector temp;
- Dune::Solvers::resizeInitializeZero(temp, v);
- f_->quadraticPart().umv(v, temp);
- temp *= 0.5;
- temp -= f_->linearPart();
- return temp*v;
- }
-
- /**
- * \brief Do one TNNMG step
- */
- void iterate() override
- {
-
- const auto& f = *f_;
- const auto& ignore = (*this->ignoreNodes_);
- auto& x = *getIterate();
-
- //std::cout << "TNNMGStep::iterate " << std::endl;
-
- //print(f.quadraticPart(), "f.quadraticPart():");
-
- //print(f.linearPart(), "f.linearPart():");
-
- //std::cout << "-- energy before smoothing: " << f(x) << std::endl;
-
- //print(x, "TNNMG iterate after smoothing:");
-
- // Nonlinear presmoothing
- for (std::size_t i=0; i<preSmoothingSteps_; ++i)
- nonlinearSmoother_->iterate();
-
- //std::cout << "- nonlinear presmoothing: success" << std::endl;
- //print(x, "TNNMG iterate after smoothing:");
-
- //std::cout << "-- energy after presmoothing: " << f(x) << std::endl;
-
- // Compute constraint/truncated linearization
- if (not linearization_)
- linearization_ = std::make_shared<Linearization>(f, ignore);
-
- linearization_->bind(x);
-
- auto&& A = linearization_->hessian();
- auto&& r = linearization_->negativeGradient();
-
- /*print(A, "TNNMG Linear problem matrix:");
- print(r, "TNNMG Linear problem rhs:");
- print(ignore, "ignore:");
- print(linearization_->truncated(), "truncation:");*/
-
- // Compute inexact solution of the linearized problem
- Dune::Solvers::resizeInitializeZero(correction_, x);
- Dune::Solvers::resizeInitializeZero(constrainedCorrection_, r);
-
-
-
- //print(constrainedCorrection_compare, "direct solver solution:");
-
- //std::cout << "- computing linear correction..." << std::endl;
-
- //linearCorrection_(A, constrainedCorrection_, r, ignore);
-
- //DUNE_THROW(Dune::Exception, "TNNMGStep: Just need to terminate here!");
-
- linearCorrection_(A, constrainedCorrection_, r, ignore);
-
- //std::cout << "- linear correction: success" << std::endl;
-
- //print(constrainedCorrection_, "contrained correction:");
-
- linearization_->extendCorrection(constrainedCorrection_, correction_);
-
- /*Vector h = x;
- h += correction_;
- std::cout << "-- energy after linear correction: " << energy(h) << std::endl;*/
-
- //std::cout << "- extended correction: success" << std::endl;
- //print(correction_, "correction:");
-
- //std::cout << "-- energy before projection: " << f(x) << std::endl;
-
- // Project onto admissible set
- projection_(f, x, correction_);
-
- //std::cout << "-- energy after projection: " << f(x) << std::endl;
-
- //std::cout << "- projection onto admissible set: success" << std::endl;
-
- //print(correction_, "correction:");
-
- double const correctionNorm = correction_.two_norm();
- correction_ /= correctionNorm;
-
-
- // Line search
- auto fv = directionalRestriction(f, x, correction_);
-
- /* std::cout << fv.quadraticPart() << " f.quadraticPart():" << std::endl;
- std::cout << fv.linearPart() << " f.linearPart():" << std::endl;
- std::cout << fv.subDifferential(0) << " f.subDifferential(0):" << std::endl;
- std::cout << fv.domain()[0] << " " << fv.domain()[0] << " f.domain():" << std::endl;
-
- std::cout << "- setup directional restriction: success" << std::endl; */
-
- dampingFactor_ = 0;
- lineSolver_(dampingFactor_, fv, false);
-
- //std::cout << "damping factor: " << dampingFactor_ << std::endl;
-
- // std::cout << "- computing damping factor: success" << std::endl;
- if (std::isnan(dampingFactor_))
- dampingFactor_ = 0;
- correction_ *= dampingFactor_;
-
- //std::cout << "damping factor: " << dampingFactor_ << std::endl;
-
- x += correction_;
-
- //std::cout << "-- energy after correction: " << energy(x) << std::endl;
- }
-
- /**
- * \brief Export the last computed damping factor
- */
- double lastDampingFactor() const
- {
- return dampingFactor_;
- }
-
- /**
- * \brief Export the last used linearization
- */
- const Linearization& linearization() const
- {
- return *linearization_;
- }
-
-private:
-
- const Functional* f_;
-
- std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother_;
- std::size_t preSmoothingSteps_ = 1;
-
- std::shared_ptr<Linearization> linearization_;
-
- //! \brief linear correction
- LinearCorrection<ConstrainedMatrix, ConstrainedVector> linearCorrection_;
-
- typename Linearization::ConstrainedVector constrainedCorrection_;
- Vector correction_;
- DefectProjection projection_;
- LineSearchSolver lineSolver_;
- double dampingFactor_;
-};
-
-
-#endif