diff --git a/dune/contact/assemblers/largedeformationcontactassembler.cc b/dune/contact/assemblers/largedeformationcontactassembler.cc
index a1bb37d63d97e43dc036e270a9765ad6669f2911..4dd8829d55bd01989a614e7a423b5df9e78c7366 100644
--- a/dune/contact/assemblers/largedeformationcontactassembler.cc
+++ b/dune/contact/assemblers/largedeformationcontactassembler.cc
@@ -7,6 +7,8 @@
#include <dune/istl/umfpack.hh>
#include <dune/matrix-vector/transpose.hh>
#include <dune/matrix-vector/addtodiagonal.hh>
+#include <dune/matrix-vector/blockmatrixview.hh>
+#include <dune/matrix-vector/algorithm.hh>
namespace Dune {
namespace Contact {
@@ -82,6 +84,8 @@ template <class GridType, class VectorType>
void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobian(const std::vector<const MatrixType*>& submat,
MatrixType& totalMatrix) const
{
+ namespace MV = Dune::MatrixVector;
+
// TODO Clean this up
// Create a block view of the grand matrix
Dune::MatrixVector::BlockMatrixView<MatrixType> bV(submat);
@@ -114,30 +118,23 @@ void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobi
}
// use fact that nonmortar and mortar matrix are dense
- auto&& it = BT_[first].begin();
- auto&& end = BT_[first].end();
-
- int offset0 = bV.col(idx0,0);
- for (; it != end; it++) {
-
+ auto offset0 = bV.col(idx0,0);
+ MV::sparseRangeFor(BT_[first], [&](auto&&, auto&& idx) {
// collected all mortar indices and nonmortar indices that are not
// part of the reduced boundary
- if (((int)it.index()>=offset0+grids_[idx0]->size(dim)) or
- ((int)it.index()<offset0) or (!nmPatch.containsVertex(it.index()-offset0)))
- contactIdx[i][1].insert(it.index());
+ if ((idx >= offset0 + this->grids_[idx0]->size(dim)) or
+ (idx < offset0) or (!nmPatch.containsVertex(idx-offset0)))
+ contactIdx[i][1].insert(idx);
else // and all nonmortar indices
- contactIdx[i][0].insert(it.index());
- }
+ contactIdx[i][0].insert(idx);
+ });
// now find all interior nonmortar nodes
- for (const auto& k : contactIdx[i][0]) {
- auto&& it = (*submat[idx0])[k].begin();
- auto&& end = (*submat[idx0])[k].end();
-
- for (; it!=end; it++)
- if (!nmPatch.containsVertex(it.index()))
- interiorIdx[i].insert(bV.row(idx0,it.index()));
- }
+ for (const auto& k : contactIdx[i][0])
+ MV::sparseRangeFor((*submat[idx0])[k], [&](auto&&, auto&& idx) {
+ if (!nmPatch.containsVertex(idx))
+ interiorIdx[i].insert(bV.row(idx0,idx));
+ });
for (size_t j=0; j<2; j++) {
@@ -167,8 +164,6 @@ void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobi
// Multiply transformation matrix to the global stiffness matrix
// ////////////////////////////////////////////////////////////////////
- namespace MV = Dune::MatrixVector;
-
indicesBABT.exportIdx(totalMatrix);
totalMatrix = 0;
@@ -188,16 +183,10 @@ void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobi
colAD[count] = 0;
- auto&& it = (*submat[idx0])[cIdx-offset0].begin();
- auto&& end = (*submat[idx0])[cIdx-offset0].end();
-
- for (; it != end; it++) {
-
- if (!nmPatch.containsVertex(it.index()))
- continue;
-
- MV::addProduct(colAD[count], *it, BT_[bV.row(idx0, it.index())][nmCol]);
- }
+ MV::sparseRangeFor((*submat[idx0])[cIdx-offset0], [&](auto&& col, auto&& idx) {
+ if (nmPatch.containsVertex(idx))
+ MV::addProduct(colAD[count], col, this->BT_[bV.row(idx0, idx)][nmCol]);
+ });
count++;
}
@@ -240,16 +229,9 @@ void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobi
colAD[count] = 0;
- auto&& it = (*submat[idx0])[cIdx-offset0].begin();
- auto&& end = (*submat[idx0])[cIdx-offset0].end();
-
- for (; it != end; it++) {
-
- if (!nmPatch.containsVertex(it.index()))
- continue;
-
- MV::addProduct(colAD[count], *it, BT_[bV.row(idx0, it.index())][mCol]);
- }
+ MV::sparseRangeFor((*submat[idx0])[cIdx-offset0], [&](auto&& col, auto&& idx) {
+ if (nmPatch.containsVertex(idx))
+ MV::addProduct(colAD[count], col, this->BT_[bV.row(idx0, idx)][mCol]);});
count++;
}
@@ -289,16 +271,11 @@ void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobi
for (size_t j=0; j<2; j++)
for (const auto& col : contactIdx[i][j]) {
- auto&& it = rowIN.begin();
- auto&& end = rowIN.end();
-
MatrixBlock sum(0);
- for (; it != end; it++) {
- if (!nmPatch.containsVertex(it.index()))
- continue;
-
- MV::addProduct(sum, *it, BT_[bV.col(idx0,it.index())][col]);
- }
+ MV::sparseRangeFor(rowIN, [&](auto&& entry, auto&& idx) {
+ if (nmPatch.containsVertex(idx))
+ MV::addProduct(sum, entry, this->BT_[bV.col(idx0, idx)][col]);
+ });
totalMatrix[intIdx][col] += sum;
MatrixBlock sumT;
@@ -327,18 +304,13 @@ void LargeDeformationContactAssembler<GridType, VectorType>::assembleExactJacobi
for(size_t j=0; j<submat[i]->N(); j++) {
// skip nonmortar rows
- if (contains(j))
- continue;
-
- const auto& mat = (*submat[i]);
-
- auto&& end = mat[j].end();
- for (auto&& it = mat[j].begin(); it != end; it++)
- if (!contains(it.index()))
- totalMatrix[bV.row(i,j)][bV.col(i,it.index())] += *it;
+ if (not contains(j))
+ MV::sparseRangeFor((*submat[i])[j], [&](auto&& col, auto&& idx) {
+ if (not contains(idx))
+ totalMatrix[bV.row(i,j)][bV.col(i,idx)] += col;
+ });
}
}
-
std::cout<<"Transformed matrix in "<<timer.stop()<<std::endl;
}
@@ -585,14 +557,11 @@ computeExactTransformationMatrix()
// finally rotate the nonmortar blocks
//TODO might not work for more than one coupling
- for(size_t rowIdx=0; rowIdx<nmNormalMatrix.N(); rowIdx++) {
+ for(size_t rowIdx=0; rowIdx<nmNormalMatrix.N(); ++rowIdx) {
int globalRowIdx = offsets[grid0Idx] + nonmortarToGlobal[i][rowIdx];
- auto& row = BT_[globalRowIdx];
-
- auto&& cEndIt = row.end();
- for(auto&& cIt = row.begin(); cIt!=cEndIt; ++cIt)
- cIt->leftmultiply(this->localCoordSystems_[globalRowIdx]);
+ for (auto& cIt : BT_[globalRowIdx])
+ cIt.leftmultiply(this->localCoordSystems_[globalRowIdx]);
}
}
}
@@ -601,6 +570,7 @@ template <class GridType, class VectorType>
void LargeDeformationContactAssembler<GridType, VectorType>::
computeLumppedTransformationMatrix()
{
+ namespace MV = Dune::MatrixVector;
/**
* Compute the mortar transformation that decouples the non-penetration constraints
* in the case of large deformations.
@@ -669,20 +639,11 @@ computeLumppedTransformationMatrix()
auto globalRowIdx = offsets[grid0Idx] + nonmortarToGlobal[i][j];
indicesBT.add(globalRowIdx,globalRowIdx);
- {
- const auto& row = nmMatrix[j];
- auto&& end = row.end();
-
- for (auto&& it = row.begin(); it != end; ++it)
- indicesBT.add(globalRowIdx, offsets[grid0Idx] + it.index());
- }
+ MV::sparseRangeFor(nmMatrix[j],[&](auto&&, auto&& idx) {
+ indicesBT.add(globalRowIdx, offsets[grid0Idx] + idx);});
- {
- const auto& row = mMatrix[j];
- auto&& end = row.end();
- for (auto&& it = row.begin(); it != end; ++it)
- indicesBT.add(globalRowIdx, offsets[grid1Idx] + it.index());
- }
+ MV::sparseRangeFor(mMatrix[j],[&](auto&&, auto&& idx) {
+ indicesBT.add(globalRowIdx, offsets[grid1Idx] + idx);});
}
}
@@ -712,27 +673,24 @@ computeLumppedTransformationMatrix()
// setup non-mortar in normal coordinates
auto normalMat = nmMat;
- for (size_t j=0; j<normalMat.N(); j++) {
-
- auto&& cEnd = normalMat[j].end();
- for(auto&& col = normalMat[j].begin(); col != cEnd; col++)
- if (nonmortarBoundary.containsVertex(col.index())) {
+ for (size_t j=0; j<normalMat.N(); j++)
+ MV::sparseRangeFor(normalMat[j], [&](auto&& col, auto&& idx) {
+ if (nonmortarBoundary.containsVertex(idx)) {
- const auto& householder = this->localCoordSystems_[offsets[grid0Idx] + col.index()];
- auto entry = (*col)[0];
- householder.mtv(entry, (*col)[0]);
+ const auto& householder = this->localCoordSystems_[offsets[grid0Idx] + idx];
+ auto entry = col[0];
+ householder.mtv(entry, col[0]);
}
- }
+ });
for (size_t j=0; j<normalMat.N(); j++) {
// invert normal component
field_type lumpedEntry(0);
- auto end = normalMat[j].end();
- for (auto it = normalMat[j].begin(); it != end; ++it)
- if (nonmortarBoundary.containsVertex(it.index()))
- lumpedEntry += (*it)[0][0];
+ MV::sparseRangeFor(normalMat[j], [&](auto&& col, auto&& idx){
+ if (nonmortarBoundary.containsVertex(idx))
+ lumpedEntry += col[0][0];});
if (std::fabs(lumpedEntry)<1e-15)
DUNE_THROW(Dune::Exception,"Warning: Normal entry of the diagonal is zero!");
@@ -744,23 +702,17 @@ computeLumppedTransformationMatrix()
// non-mortar entries
auto globalRowIdx = offsets[grid0Idx] + nonmortarToGlobal[i][j];
- {
- auto end = normalMat[j].end();
- for (auto it = normalMat[j].begin(); it != end; ++it) {
- BT_[globalRowIdx][offsets[grid0Idx] + it.index()][0].axpy(lumpedEntry,(*it)[0]);
- if (nonmortarBoundary.containsVertex(it.index()))
- BT_[globalRowIdx][offsets[grid0Idx] + it.index()][0][0] = 0;
- }
+ MV::sparseRangeFor(normalMat[j], [&](auto&& col, auto&& idx) {
+ this->BT_[globalRowIdx][offsets[grid0Idx] + idx][0].axpy(lumpedEntry,col[0]);
+ if (nonmortarBoundary.containsVertex(idx))
+ this->BT_[globalRowIdx][offsets[grid0Idx] + idx][0][0] = 0;
+ });
- BT_[globalRowIdx][globalRowIdx][0][0] = lumpedEntry;
- }
+ BT_[globalRowIdx][globalRowIdx][0][0] = lumpedEntry;
- {
- const auto& row = mMat[j];
- auto&& end = row.end();
- for (auto&& it = row.begin(); it != end; ++it)
- BT_[globalRowIdx][offsets[grid1Idx] + it.index()][0].axpy(lumpedEntry,(*it)[0]);
- }
+ MV::sparseRangeFor(mMat[j], [&](auto&& col, auto&& idx) {
+ this->BT_[globalRowIdx][offsets[grid1Idx] + idx][0].axpy(lumpedEntry, col[0]);
+ });
}
// finally rotate the nonmortar blocks
@@ -780,6 +732,7 @@ template <class GridType, class VectorType>
void LargeDeformationContactAssembler<GridType, VectorType>::
computeInverseTransformationMatrix()
{
+ namespace MV = Dune::MatrixVector;
/**
* Compute the inverse of the mortar transformation that decouples the non-penetration constraints
* in the case of large deformations.
@@ -890,60 +843,46 @@ computeInverseTransformationMatrix()
// non-mortar entries
auto globalRowIdx = offsets[grid0Idx] + nonmortarToGlobal[i][j];
- {
- const auto& row = nmMat[j];
- auto&& end = row.end();
-
- // TODO: only add entries that are non-zero in normal coords
- for (auto&& it = row.begin(); it != end; ++it) {
-
- MatrixBlock invTrans(0);
-
- // for non-mortar components, the entrie is twice multiplied by the Householder
- // reflection, hence we don't have to compute normal coordinates
- if (nonmortarBoundary.containsVertex(it.index()) and
- nonmortarToGlobal[i][j]== it.index())
- invTrans = this->localCoordSystems_[offsets[grid0Idx] + it.index()];
+ // TODO: only add entries that are non-zero in normal coords
+ MV::sparseRangeFor(nmMat[j], [&](auto&& col, auto&& idx) {
+ MatrixBlock invTrans(0);
- invTrans[0] = (*it)[0];
- invTrans[0] *= -1;
+ // for non-mortar components, the entrie is twice multiplied by the Householder
+ // reflection, hence we don't have to compute normal coordinates
+ if (nonmortarBoundary.containsVertex(idx) and
+ nonmortarToGlobal[i][j]== idx)
+ invTrans = this->localCoordSystems_[offsets[grid0Idx] + idx];
- inverse_[globalRowIdx][offsets[grid0Idx]+it.index()] = invTrans;
- }
- }
+ invTrans[0] = col[0];
+ invTrans[0] *= -1;
- {
- const auto& row = mMat[j];
- auto&& end = row.end();
- for (auto&& it = row.begin(); it != end; ++it) {
+ inverse_[globalRowIdx][offsets[grid0Idx] + idx] = invTrans;
+ });
- MatrixBlock invTrans(0);
- invTrans[0] -= (*it)[0];
- inverse_[offsets[grid0Idx]+nonmortarToGlobal[i][j]][offsets[grid1Idx]+it.index()] = invTrans;
- }
- }
+ MV::sparseRangeFor(mMat[j], [&](auto&& col, auto&& idx) {
+ MatrixBlock invTrans(0);
+ invTrans[0] -= col[0];
+ inverse_[globalRowIdx][offsets[grid1Idx] + idx] = invTrans;
+ });
}
}
// compute the norms of the columns, these determine the length of the
// transformed basis vectors and are taken into account by the trust-region
- const auto& inverse = inverse_;
- colNorms_.resize(inverse.N());
+ colNorms_.resize(inverse_.N());
colNorms_ = 0;
- for (size_t i=0; i<inverse.N(); i++) {
- auto cEnd = inverse[i].end();
- for (auto col = inverse[i].begin(); col != cEnd; col++)
- for (int j=0; j<dim;j++)
- for (int k=0; k<dim;k++)
- colNorms_[col.index()][k] += std::pow((*col)[j][k],2);
- }
-
- field_type min(std::numeric_limits<field_type>::max()),max(-1);
- for (size_t i=0; i<inverse.N(); i++)
- for (int j=0; j<dim;j++) {
- field_type root = std::sqrt(colNorms_[i][j]);
- colNorms_[i][j] = 1.0/root;
+ MV::sparseMatrixFor(inverse_, [&](auto&& entry, auto&&, auto&& colIdx) {
+ for (int j=0; j < dim; ++j)
+ for (int k=0; k < dim; ++k)
+ colNorms_[colIdx][k] += std::pow(entry[j][k], 2);
+ });
+
+ field_type min(std::numeric_limits<field_type>::max()), max(-1);
+ for (auto& colNorm : colNorms_)
+ for (int j=0; j < dim; ++j) {
+ field_type root = std::sqrt(colNorm[j]);
+ colNorm[j] = 1.0/root;
min = std::fmin(min, root);
max = std::fmax(max, root);
@@ -951,13 +890,14 @@ computeInverseTransformationMatrix()
std::cout<<" Transformed basis vector norms range from "<<min<<" to "<<max<<std::endl;
}
-
template <class GridType, class VectorType>
template <class VectorTypeContainer>
void LargeDeformationContactAssembler<GridType, VectorType>::
assembleExactRightHandSide(const VectorTypeContainer& rhs,
VectorType& totalRhs) const
{
+ namespace MV = Dune::MatrixVector;
+
// first combine all regular parts
this->concatenateVectors(rhs,totalRhs);
@@ -991,14 +931,10 @@ assembleExactRightHandSide(const VectorTypeContainer& rhs,
// indices of normal nm matrix
Dune::MatrixIndexSet nmIsN(nmMatrix.N(),nmMatrix.N());
// for the normal projection the tangential matrix is not square
- for (size_t j=0; j<nmMatrix.N(); j++) {
- const auto& row = nmMatrix[j];
- auto&& end = row.end();
-
- for (auto&& it = row.begin(); it != end; ++it)
- if (nonmortarBoundary.containsVertex(it.index()))
- nmIsN.add(globalToLocal[it.index()],j);
- }
+ MV::sparseMatrixFor(nmMatrix, [&](auto&&, auto&& rowIdx, auto&& colIdx) {
+ if (nonmortarBoundary.containsVertex(colIdx))
+ nmIsN.add(globalToLocal[colIdx], rowIdx);
+ });
using ScalarMatrix = BCRSMatrix<FieldMatrix<field_type,1 , 1> >;
ScalarMatrix nmMatN;
@@ -1028,22 +964,17 @@ assembleExactRightHandSide(const VectorTypeContainer& rhs,
totalRhs[globalRowIdx][0] = 0;
}
- for (size_t j=0; j<nmMatT.N(); j++) {
- auto& row = nmMatT[j];
- for(auto&& col = row.begin(); col != nmMatT[j].end(); col++) {
+ MV::sparseMatrixFor(nmMatT, [&](auto&& entry, auto&& rowIdx, auto&& colIdx){
+ if (nonmortarBoundary.containsVertex(colIdx)) {
+ const auto& householder = this->localCoordSystems_[offsets[gridIdx[0]] + colIdx];
- if (!nonmortarBoundary.containsVertex(col.index()))
- continue;
- const auto& householder = this->localCoordSystems_[offsets[gridIdx[0]] + col.index()];
-
- auto entry = (*col)[0];
- householder.mtv(entry, (*col)[0]);
+ auto val = entry[0];
+ householder.mtv(val, entry[0]);
- // if (std::fabs((*col)[0][0])>1e-15)
- nmMatN[globalToLocal[col.index()]][j][0][0] = (*col)[0][0];
- (*col)[0][0] = 0;
+ nmMatN[globalToLocal[colIdx]][rowIdx][0][0] = entry[0][0];
+ entry[0][0] = 0;
}
- }
+ });
// solve the small system nmMatN*transRhsNm = rhsNm;
auto transRhsNm = rhsNm;
@@ -1060,16 +991,12 @@ assembleExactRightHandSide(const VectorTypeContainer& rhs,
std::vector<std::reference_wrapper<const JacobianType> > linMatrices;
linMatrices = {{std::cref(nmMatT), std::cref(coupling->exactMortarMatrix())}};
- for (size_t j =0; j<2; j++) {
- const auto& mat = linMatrices[j].get();
+ for (size_t j =0; j < 2; ++j)
+ MV::sparseMatrixFor(linMatrices[j].get(), [&](auto&& entry, auto&& rowIdx, auto&& colIdx) {
+ totalRhs[offsets[gridIdx[j]] + colIdx].axpy(transRhsNm[rowIdx][0], entry[0]);});
- for (size_t p=0; p<linMatrices[j].get().N(); p++)
- for (auto col = mat[p].begin(); col != mat[p].end();col++)
- totalRhs[offsets[gridIdx[j]]+col.index()].axpy(transRhsNm[p][0],(*col)[0]);
- }
-
- for (size_t p=0; p<transRhsNm.size(); p++)
- totalRhs[offsets[gridIdx[0]]+nonmortarToGlobal[p]][0] += transRhsNm[p][0];
+ for (size_t p=0; p < transRhsNm.size(); ++p)
+ totalRhs[offsets[gridIdx[0]] + nonmortarToGlobal[p]][0] += transRhsNm[p][0];
}
}
@@ -1078,6 +1005,8 @@ template <class VectorTypeContainer>
void LargeDeformationContactAssembler<GridType, VectorType>::
postprocess(const VectorType& totalX, VectorTypeContainer& x) const
{
+ namespace MV = Dune::MatrixVector;
+
auto untransformedX = totalX;
// compute offsets for the involved grids
@@ -1110,14 +1039,10 @@ postprocess(const VectorType& totalX, VectorTypeContainer& x) const
// indices of normal nm matrix
Dune::MatrixIndexSet nmIsN(nmMatrix.N(),nmMatrix.N());
// for the normal projection the tangential matrix is not square
- for (size_t j=0; j<nmMatrix.N(); j++) {
- const auto& row = nmMatrix[j];
- auto&& end = row.end();
-
- for (auto&& it = row.begin(); it != end; ++it)
- if (nonmortarBoundary.containsVertex(it.index()))
- nmIsN.add(j, globalToLocal[it.index()]);
- }
+ MV::sparseMatrixFor(nmMatrix, [&](auto&&, auto&& rowIdx, auto&& colIdx) {
+ if (nonmortarBoundary.containsVertex(colIdx))
+ nmIsN.add(rowIdx, globalToLocal[colIdx]);
+ });
using ScalarMatrix = BCRSMatrix<FieldMatrix<field_type,1 , 1> >;
ScalarMatrix nmMatN;
@@ -1156,9 +1081,8 @@ postprocess(const VectorType& totalX, VectorTypeContainer& x) const
// add mortar terms
const auto& mMatrix = coupling->exactMortarMatrix();
- for (size_t j=0; j<mMatrix.N(); j++)
- for(auto col = mMatrix[j].begin(); col != mMatrix[j].end(); col++)
- rhsNm[j] += (*col)[0]*totalX[offsets[gridIdx[1]] + col.index()];
+ MV::sparseMatrixFor(mMatrix, [&](auto&& entry, auto&& rowIdx, auto&& colIdx) {
+ rhsNm[rowIdx] += entry[0]*totalX[offsets[gridIdx[1]] + colIdx];});
// we actually need the negative term
rhsNm *= -1;