From 9a075dfedb4d25f2b8d11271d3119efd9cdec573 Mon Sep 17 00:00:00 2001
From: Jonathan Youett <youett@math.fu-berlin.de>
Date: Tue, 23 Feb 2016 15:25:16 +0100
Subject: [PATCH] Mandel restriction now also works if DenseMultigridTransfer
is used.
By introducing a template helper method that gets a possibly blocked
restriction matrix, this class can now also be used together with the
DenseMultigridStep. This is used for example in the case of large
deformation contact problems.
---
.../transferoperators/mandelobsrestrictor.cc | 47 +++++++++++++------
.../transferoperators/mandelobsrestrictor.hh | 12 +++++
2 files changed, 44 insertions(+), 15 deletions(-)
diff --git a/dune/solvers/transferoperators/mandelobsrestrictor.cc b/dune/solvers/transferoperators/mandelobsrestrictor.cc
index 2593b433..50f31679 100644
--- a/dune/solvers/transferoperators/mandelobsrestrictor.cc
+++ b/dune/solvers/transferoperators/mandelobsrestrictor.cc
@@ -2,6 +2,7 @@
// vi: set et ts=8 sw=4 sts=4:
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+#include <dune/solvers/transferoperators/densemultigridtransfer.hh>
template <class VectorType>
void MandelObstacleRestrictor<VectorType>::
@@ -12,43 +13,59 @@ restrict(const std::vector<BoxConstraint<field_type, blocksize> >&f,
const MultigridTransfer<VectorType>& transfer,
const Dune::BitSetVector<blocksize>& critical)
{
- if (!(dynamic_cast<const CompressedMultigridTransfer<VectorType>*>(&transfer)))
- DUNE_THROW(Dune::NotImplemented, "Implementation assumes transfer to be of type CompressedMultigridTransfer!");
- const Dune::BCRSMatrix<Dune::FieldMatrix<field_type, 1, 1> >& transferMat =
- dynamic_cast<const CompressedMultigridTransfer<VectorType>*>(&transfer)->getMatrix();
+ if (dynamic_cast<const CompressedMultigridTransfer<VectorType>*>(&transfer))
+ restrict(f,t,fHasObstacle,tHasObstacle,
+ dynamic_cast<const CompressedMultigridTransfer<VectorType>*>(&transfer)->getMatrix(),critical);
+ else if (dynamic_cast<const DenseMultigridTransfer<VectorType>*>(&transfer))
+ restrict(f,t,fHasObstacle,tHasObstacle,
+ dynamic_cast<const DenseMultigridTransfer<VectorType>*>(&transfer)->getMatrix(),critical);
+ else
+ DUNE_THROW(Dune::NotImplemented, "Implementation assumes transfer to be of type CompressedMultigridTransfer"
+ "or DenseMultigridTransfer!");
+}
+
+template <class VectorType>
+template <int dim>
+void MandelObstacleRestrictor<VectorType>::
+restrict(const std::vector<BoxConstraint<field_type, blocksize> >&f,
+ std::vector<BoxConstraint<field_type, blocksize> > &t,
+ const Dune::BitSetVector<blocksize>& fHasObstacle,
+ const Dune::BitSetVector<blocksize>& tHasObstacle,
+ const Dune::BCRSMatrix<Dune::FieldMatrix<field_type,dim,dim> >& transferMat,
+ const Dune::BitSetVector<blocksize>& critical)
+{
+ // dim should be either 1 for the compressed transfer or dim
+ assert(dim==blocksize || dim==1);
t.resize(transferMat.M());
for (size_t i=0; i<t.size(); i++)
t[i].clear();
// We use the multigrid transfer matrix to find all relevant fine grid dofs
- typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> >::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
// Loop over all coarse grid dofs
for (size_t i=0; i<transferMat.N(); i++) {
if (fHasObstacle[i].none())
continue;
- const RowType& row = transferMat[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
+ const auto& row = transferMat[i];
// Loop over all fine grid base functions in the support of tObsIt
- for(; cIt!=cEndIt; ++cIt) {
+ for(auto cIt=row.begin(); cIt != row.end(); ++cIt) {
if (tHasObstacle[cIt.index()].none())
continue;
// Sort out spurious entries due to numerical dirt
- if ((*cIt)[0][0] < 1e-3)
- continue;
+ // initialise with the first entry to ensure to handle the case dim==1
+ std::vector<bool> truncated(blocksize,((*cIt)[0][0]<1e-3));
+ for (int j=1; j<dim; j++)
+ truncated[j] = ((*cIt)[j][j] < 1e-3);
+
for (int j=0; j<blocksize; j++)
- if (!critical[i][j] && fHasObstacle[i][j] && tHasObstacle[cIt.index()][j]) {
+ if (not truncated[j] and not critical[i][j]) {
t[cIt.index()].lower(j) = std::max(t[cIt.index()].lower(j), f[i].lower(j));
t[cIt.index()].upper(j) = std::min(t[cIt.index()].upper(j), f[i].upper(j));
}
diff --git a/dune/solvers/transferoperators/mandelobsrestrictor.hh b/dune/solvers/transferoperators/mandelobsrestrictor.hh
index 5a41c098..03f0f8dc 100644
--- a/dune/solvers/transferoperators/mandelobsrestrictor.hh
+++ b/dune/solvers/transferoperators/mandelobsrestrictor.hh
@@ -4,6 +4,8 @@
#define DUNE_SOLVERS_TRANSFEROPERATORS_MANDEL_OBSTACLE_RESTRICTOR_HH
#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
#include <dune/solvers/transferoperators/obstaclerestrictor.hh>
#include <dune/solvers/transferoperators/multigridtransfer.hh>
#include <dune/solvers/common/boxconstraint.hh>
@@ -32,6 +34,16 @@ public:
const Dune::BitSetVector<blocksize>& tHasObstacle,
const MultigridTransfer<VectorType>& transfer,
const Dune::BitSetVector<blocksize>& critical);
+
+private:
+ //! Helper method that allows to use the restrictor for both compressed and cense transfer operator
+ template <int dim>
+ void restrict(const std::vector<BoxConstraint<field_type,blocksize> >& f,
+ std::vector<BoxConstraint<field_type,blocksize> >& t,
+ const Dune::BitSetVector<blocksize>& fHasObstacle,
+ const Dune::BitSetVector<blocksize>& tHasObstacle,
+ const Dune::BCRSMatrix<Dune::FieldMatrix<field_type,dim,dim> >& transferMat,
+ const Dune::BitSetVector<blocksize>& critical);
};
// include implementation
--
GitLab