From f3abb6005befbb853ca8be631227d642a2d43fc9 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Wed, 8 May 2019 19:04:29 +0200
Subject: [PATCH] .
---
dune/tectonic/globalfriction.hh | 2 +-
dune/tectonic/globalfrictiondata.hh | 1 +
dune/tectonic/globalratestatefriction.hh | 125 ++++--
.../transformedglobalratestatefriction.hh | 129 ------
src/CMakeLists.txt | 6 +-
src/assemblers.cc | 43 +-
src/assemblers.hh | 7 +
src/assemblers_tmpl.cc | 16 +-
src/boundarycondition.hh | 42 --
src/data-structures/body.cc | 69 +---
src/data-structures/body.hh | 24 +-
src/data-structures/levelcontactnetwork.cc | 120 +++++-
src/data-structures/levelcontactnetwork.hh | 34 +-
.../levelcontactnetwork_tmpl.cc | 2 +
src/data-structures/program_state.hh | 2 +-
src/factories/levelcontactnetworkfactory.hh | 4 +-
src/factories/stackedblocksfactory.cc | 25 +-
src/frictioncouplingpair.hh | 37 ++
src/io/hdf5/frictionalboundary-writer.cc | 2 +-
src/multi-body-problem-data/weakpatch.hh | 2 +-
src/multi-body-problem.cc | 25 +-
src/spatial-solving/fixedpointiterator.cc | 255 ++----------
src/spatial-solving/fixedpointiterator.hh | 19 +-
.../fixedpointiterator_tmpl.cc | 31 +-
src/spatial-solving/solverfactory.hh | 2 +-
src/spatial-solving/solverfactory_tmpl.cc | 11 +-
src/tests/contactmerge.cc | 40 ++
src/tests/couplingtest.hh | 192 +++++++++
src/tests/gridgluefrictiontest.cc | 317 +++++++++------
src/tests/nonoverlappingcouplingtest.cc | 375 ++++++++++++++++++
src/time-stepping/adaptivetimestepper.cc | 20 +-
src/time-stepping/adaptivetimestepper.hh | 13 +-
src/time-stepping/adaptivetimestepper_tmpl.cc | 20 +-
src/time-stepping/coupledtimestepper.cc | 7 +-
src/time-stepping/coupledtimestepper.hh | 10 +-
src/time-stepping/coupledtimestepper_tmpl.cc | 20 +-
src/time-stepping/state.cc | 35 +-
src/time-stepping/state.hh | 9 +-
.../state/ageinglawstateupdater.cc | 85 ++--
.../state/ageinglawstateupdater.hh | 38 +-
.../state/sliplawstateupdater.cc | 80 ++--
.../state/sliplawstateupdater.hh | 40 +-
src/time-stepping/state/stateupdater.hh | 66 ++-
src/time-stepping/state_tmpl.cc | 29 +-
src/utils/almostequal.hh | 2 +
src/utils/debugutils.hh | 1 +
46 files changed, 1466 insertions(+), 968 deletions(-)
delete mode 100644 dune/tectonic/transformedglobalratestatefriction.hh
create mode 100644 src/frictioncouplingpair.hh
create mode 100644 src/tests/contactmerge.cc
create mode 100644 src/tests/couplingtest.hh
create mode 100644 src/tests/nonoverlappingcouplingtest.cc
diff --git a/dune/tectonic/globalfriction.hh b/dune/tectonic/globalfriction.hh
index e87fc139..83395851 100644
--- a/dune/tectonic/globalfriction.hh
+++ b/dune/tectonic/globalfriction.hh
@@ -84,7 +84,7 @@ class GlobalFriction { //: public Dune::TNNMG::NonsmoothConvexFunctional<> {
return ret;
}
- void virtual updateAlpha(ScalarVector const &alpha) = 0;
+ void virtual updateAlpha(const std::vector<ScalarVector>& alpha) = 0;
};
#endif
diff --git a/dune/tectonic/globalfrictiondata.hh b/dune/tectonic/globalfrictiondata.hh
index f3851f99..0fffac2c 100644
--- a/dune/tectonic/globalfrictiondata.hh
+++ b/dune/tectonic/globalfrictiondata.hh
@@ -27,6 +27,7 @@ template <int dimension> class GlobalFrictionData {
Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
Dune::FieldVector<double, 1>>;
+public:
double virtual const &C() const = 0;
double virtual const &L() const = 0;
double virtual const &V0() const = 0;
diff --git a/dune/tectonic/globalratestatefriction.hh b/dune/tectonic/globalratestatefriction.hh
index e233523d..edfb05e0 100644
--- a/dune/tectonic/globalratestatefriction.hh
+++ b/dune/tectonic/globalratestatefriction.hh
@@ -10,62 +10,123 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfrictiondata.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/index-in-sorted-range.hh>
-template <class Matrix, class Vector, class ScalarFriction, class GridView>
+#include "../../src/frictioncouplingpair.hh"
+
+template <class Matrix, class Vector, class ScalarFriction, class GridType>
class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
public:
using GlobalFriction<Matrix, Vector>::block_size;
using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
private:
- using typename GlobalFriction<Matrix, Vector>::ScalarVector;
+ using Base = GlobalFriction<Matrix, Vector>;
+
+ using field_type = typename Vector::field_type;
+ using typename Base::ScalarVector;
+ using typename Base::LocalVectorType;
+
+ using FrictionCouplingPair = FrictionCouplingPair<GridType, LocalVectorType, field_type>;
+ using ContactCoupling = Dune::Contact::DualMortarCoupling<field_type, GridType>;
+
+ size_t bodyIndex(const size_t globalIdx) {
+ size_t i=0;
+
+ for (; i<maxIndex_.size(); i++) {
+ if (globalIdx < maxIndex_[i])
+ break;
+ }
+ return i;
+ }
public:
- GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
- GlobalFrictionData<block_size> const &frictionInfo,
- ScalarVector const &weights,
- ScalarVector const &weightedNormalStress)
- : restrictions(), localToGlobal(), zeroFriction() {
- auto const gridView = frictionalBoundary.gridView();
- Dune::MultipleCodimMultipleGeomTypeMapper<
- GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
- for (auto it = gridView.template begin<block_size>();
- it != gridView.template end<block_size>(); ++it) {
- auto const i = vertexMapper.index(*it);
-
- if (not frictionalBoundary.containsVertex(i))
- continue;
-
- localToGlobal.emplace_back(i);
- restrictions.emplace_back(weights[i], weightedNormalStress[i],
- frictionInfo(geoToPoint(it->geometry())));
+ GlobalRateStateFriction(const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+ const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings, // contains frictionInfo
+ const std::vector<ScalarVector>& weights,
+ const std::vector<ScalarVector>& weightedNormalStress)
+ : restrictions_(), localToGlobal_(), zeroFriction_() {
+
+ assert(contactCouplings.size() == couplings.size());
+ assert(weights.size() == weightedNormalStress.size());
+
+ const auto nBodies = weights.size();
+
+ std::vector<std::vector<int>> nonmortarBodies(nBodies); // first index body, second index coupling
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ const auto nonmortarIdx = couplings[i]->gridIdx_[0];
+ nonmortarBodies[nonmortarIdx].emplace_back(i);
+ }
+
+ maxIndex_.resize(nBodies);
+ size_t globalIdx = 0;
+
+ for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+ const auto& couplingIndices = nonmortarBodies[bodyIdx];
+
+ if (couplingIndices.size()==0)
+ continue;
+
+ const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
+
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+ for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
+ const auto i = vertexMapper.index(*it);
+
+ for (size_t j=0; j<couplingIndices.size(); j++) {
+ const auto couplingIdx = couplingIndices[j];
+
+ if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
+ continue;
+
+ localToGlobal_.emplace_back(i);
+ restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
+ couplings[i]->frictionData()(geoToPoint(it->geometry())));
+ break;
+ }
+
+ globalIdx++;
+ }
+ maxIndex_[bodyIdx] = globalIdx;
}
- assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
- assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
}
- void updateAlpha(ScalarVector const &alpha) override {
- for (size_t j = 0; j < restrictions.size(); ++j)
- restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
+ void updateAlpha(const std::vector<ScalarVector>& alpha) override {
+ for (size_t j = 0; j < restrictions_.size(); ++j) {
+ const auto bodyIdx = bodyIndex(j);
+ const auto globalDof = localToGlobal_[j];
+
+ size_t bodyDof;
+ if (bodyIdx>0) {
+ bodyDof = globalDof - maxIndex_[bodyIdx-1];
+ } else {
+ bodyDof = globalDof;
+ }
+ restrictions_[j].updateAlpha(alpha[bodyIdx][bodyDof]);
+ }
}
/*
Return a restriction of the outer function to the i'th node.
*/
LocalNonlinearity const &restriction(size_t i) const override {
- auto const index = indexInSortedRange(localToGlobal, i);
- if (index == localToGlobal.size())
- return zeroFriction;
- return restrictions[index];
+ auto const index = indexInSortedRange(localToGlobal_, i);
+ if (index == localToGlobal_.size())
+ return zeroFriction_;
+ return restrictions_[index];
}
private:
- std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
- std::vector<size_t> localToGlobal;
- WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
+ std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions_;
+ std::vector<size_t> maxIndex_; // max global index per body
+ std::vector<size_t> localToGlobal_;
+ WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction_;
};
#endif
diff --git a/dune/tectonic/transformedglobalratestatefriction.hh b/dune/tectonic/transformedglobalratestatefriction.hh
deleted file mode 100644
index 0a4cd4a6..00000000
--- a/dune/tectonic/transformedglobalratestatefriction.hh
+++ /dev/null
@@ -1,129 +0,0 @@
-#ifndef DUNE_TECTONIC_TRANSFORMED_GLOBALRATESTATEFRICTION_HH
-#define DUNE_TECTONIC_TRANSFORMED_GLOBALRATESTATEFRICTION_HH
-
-#include <vector>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/contact/assemblers/nbodyassembler.hh>
-
-//#include <dune/tectonic/geocoordinate.hh>
-//#include <dune/tectonic/globalfrictiondata.hh>
-#include <dune/tectonic/globalratestatefriction.hh>
-//#include <dune/tectonic/index-in-sorted-range.hh>
-
-template <class Matrix, class Vector, class ScalarFriction, class GridView>
-class TransformedGlobalRateStateFriction : public GlobalRateStateFriction<Matrix, Vector, ScalarFriction, GridView> {
-public:
- using NBodyAssembler = Dune::Contact::NBodyAssembler<typename GridView::Grid, Vector>;
-
- using GlobalFriction<Matrix, Vector>::block_size;
- using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
-
-private:
- using typename GlobalFriction<Matrix, Vector>::ScalarVector;
-
-public:
- TransformedGlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
- GlobalFrictionData<block_size> const &frictionInfo,
- ScalarVector const &weights,
- ScalarVector const &weightedNormalStress)
- : restrictions(), localToGlobal(), zeroFriction() {
-
- auto const gridView = frictionalBoundary.gridView();
- Dune::MultipleCodimMultipleGeomTypeMapper<
- GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
- for (auto it = gridView.template begin<block_size>();
- it != gridView.template end<block_size>(); ++it) {
- auto const i = vertexMapper.index(*it);
-
- if (not frictionalBoundary.containsVertex(i))
- continue;
-
- localToGlobal.emplace_back(i);
- restrictions.emplace_back(weights[i], weightedNormalStress[i],
- frictionInfo(geoToPoint(it->geometry())));
- }
-
- assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
- assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
- }
-
- double operator()(Vector const &x) const {
- Vector nodalX(x.size());
- nBodyAssembler_.
- double tmp = 0;
- for (size_t i = 0; i < x.size(); ++i) {
- tmp += restriction(i)(x[i]);
- }
- return tmp;
- }
-
- /*
- Return a restriction of the outer function to the i'th node.
- */
- LocalNonlinearity const virtual &restriction(size_t i) const = 0;
-
- void addHessian(Vector const &v, Matrix &hessian) const {
- for (size_t i = 0; i < v.size(); ++i)
- restriction(i).addHessian(v[i], hessian[i][i]);
- }
-
- void directionalDomain(Vector const &, Vector const &,
- Dune::Solvers::Interval<double> &dom) const {
- dom[0] = -std::numeric_limits<double>::max();
- dom[1] = std::numeric_limits<double>::max();
- }
-
- void directionalSubDiff(
- Vector const &u, Vector const &v,
- Dune::Solvers::Interval<double> &subdifferential) const {
- subdifferential[0] = subdifferential[1] = 0;
- for (size_t i = 0; i < u.size(); ++i) {
- Dune::Solvers::Interval<double> D;
- restriction(i).directionalSubDiff(u[i], v[i], D);
- subdifferential[0] += D[0];
- subdifferential[1] += D[1];
- }
- }
-
- void addHessianIndices(Dune::MatrixIndexSet &indices) const {
- for (size_t i = 0; i < indices.rows(); ++i)
- indices.add(i, i);
- }
-
- void addGradient(Vector const &v, Vector &gradient) const {
- for (size_t i = 0; i < v.size(); ++i)
- restriction(i).addGradient(v[i], gradient[i]);
- }
-
- double regularity(size_t i, typename Vector::block_type const &x) const {
- return restriction(i).regularity(x);
- }
-
- ScalarVector coefficientOfFriction(Vector const &x) const {
- ScalarVector ret(x.size());
- for (size_t i = 0; i < x.size(); ++i)
- ret[i] = restriction(i).coefficientOfFriction(x[i]);
- return ret;
- }
-
- /*
- Return a restriction of the outer function to the i'th node.
- */
- LocalNonlinearity const &restriction(size_t i) const override {
- auto const index = indexInSortedRange(localToGlobal, i);
- if (index == localToGlobal.size())
- return zeroFriction;
- return restrictions[index];
- }
-
-private:
- const NBodyAssembler& nBodyAssembler_;
-};
-#endif
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 32c35011..7c012be4 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -25,7 +25,7 @@ set(SW_SOURCE_FILES
set(MSW_SOURCE_FILES
assemblers.cc
- nodalweights.cc
+ #nodalweights.cc
data-structures/body.cc
#data-structures/contactnetwork.cc
data-structures/enumparser.cc
@@ -45,9 +45,9 @@ set(MSW_SOURCE_FILES
multi-body-problem-data/grid/mygrids.cc
multi-body-problem-data/grid/simplexmanager.cc
multi-body-problem.cc
- spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc
- spatial-solving/solverfactory_old.cc
+ spatial-solving/fixedpointiterator.cc
+ #spatial-solving/solverfactory_old.cc
time-stepping/adaptivetimestepper.cc
time-stepping/coupledtimestepper.cc
time-stepping/rate.cc
diff --git a/src/assemblers.cc b/src/assemblers.cc
index d19b8e32..1160fd59 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -12,7 +12,6 @@
#include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functiontools/p0p1interpolation.hh>
@@ -31,6 +30,16 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView) :
cellAssembler(cellBasis, cellBasis),
vertexAssembler(vertexBasis, vertexBasis) {}
+template <class GridView, int dimension>
+template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
+void MyAssembler<GridView, dimension>::assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
+ GlobalVectorType& b,
+ const BoundaryPatch<GridView>& boundaryPatch,
+ bool initializeVector) const {
+
+ vertexAssembler.assembleBoundaryFunctional(localAssembler, b, boundaryPatch, initializeVector);
+}
+
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
BoundaryPatch<GridView> const &frictionalBoundary,
@@ -146,38 +155,6 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
}
-template <class GridView, int dimension>
-auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
- Config::FrictionModel frictionModel,
- BoundaryPatch<GridView> const &frictionalBoundary,
- GlobalFrictionData<dimension> const &frictionInfo,
- ScalarVector const &weightedNormalStress) const
--> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
-
- // Lumping of the nonlinearity
- ScalarVector weights;
- {
- NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
- frictionalBoundaryAssembler(std::make_shared<
- ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
- 1));
- vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
- weights, frictionalBoundary);
- }
- switch (frictionModel) {
- case Config::Truncated:
- return std::make_shared<GlobalRateStateFriction<
- Matrix, Vector, TruncatedRateState, GridView>>(
- frictionalBoundary, frictionInfo, weights, weightedNormalStress);
- case Config::Regularised:
- return std::make_shared<GlobalRateStateFriction<
- Matrix, Vector, RegularisedRateState, GridView>>(
- frictionalBoundary, frictionInfo, weights, weightedNormalStress);
- default:
- assert(false);
- }
-}
-
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleVonMisesStress(
double youngModulus,
diff --git a/src/assemblers.hh b/src/assemblers.hh
index a32870a9..695a94af 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -12,6 +12,7 @@
#include <dune/fufem/functionspacebases/p0basis.hh>
#pragma clang diagnostic pop
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/boundarypatch.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/globalfrictiondata.hh>
@@ -48,6 +49,12 @@ class MyAssembler {
public:
MyAssembler(GridView const &gridView);
+ template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
+ void assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
+ GlobalVectorType& b,
+ const BoundaryPatch<GridView>& boundaryPatch,
+ bool initializeVector=true) const;
+
void assembleFrictionalBoundaryMass(
BoundaryPatch<GridView> const &frictionalBoundary,
ScalarMatrix &frictionalBoundaryMass) const;
diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index edc0a041..acf72e7a 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -5,6 +5,18 @@
#include "explicitgrid.hh"
#include "explicitvectors.hh"
-#include "data-structures/body.hh"
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
-template class MyAssembler<Body<Grid, Vector, MY_DIM>::DeformedLeafGridView, MY_DIM>;
+#include "assemblers.hh"
+
+template class MyAssembler<DefLeafGridView, MY_DIM>;
+
+
+using MyNeumannBoundaryAssembler = NeumannBoundaryAssembler<DeformedGrid, typename ScalarVector::block_type>;
+
+template void MyAssembler<DefLeafGridView, MY_DIM>::assembleBoundaryFunctional<MyNeumannBoundaryAssembler, ScalarVector>(
+ MyNeumannBoundaryAssembler& localAssembler,
+ ScalarVector& b,
+ const BoundaryPatch<DefLeafGridView>& boundaryPatch,
+ bool initializeVector) const;
diff --git a/src/boundarycondition.hh b/src/boundarycondition.hh
index 2ab5f2b7..088ebc97 100644
--- a/src/boundarycondition.hh
+++ b/src/boundarycondition.hh
@@ -71,46 +71,4 @@ class BoundaryCondition {
}
};
-
-#include <dune/fufem/geometry/convexpolyhedron.hh>
-
-#include <dune/tectonic/globalfrictiondata.hh>
-
-template <class GridView, class LocalVectorType, int dims>
-class FrictionBoundaryCondition : public BoundaryCondition<GridView, dims>{
-private:
- using Base = BoundaryCondition<GridView, dims>;
- using LocalVector = LocalVectorType;
-
- // friction data
- std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
- std::shared_ptr<GlobalFrictionData<dims>> frictionData_;
-
-public:
- FrictionBoundaryCondition() :
- Base("friction")
- {}
-
- FrictionBoundaryCondition(std::shared_ptr<typename Base::BoundaryPatch> patch, std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch, std::shared_ptr<GlobalFrictionData<dims>> frictionData) :
- Base(patch, nullptr, "friction"),
- weakeningPatch_(weakeningPatch),
- frictionData_(frictionData)
- {}
-
- void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
- weakeningPatch_ = weakeningPatch;
- }
-
- void setFrictionData(std::shared_ptr<GlobalFrictionData<dims>> frictionData) {
- frictionData_ = frictionData;
- }
-
- const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch() const {
- return weakeningPatch_;
- }
-
- const std::shared_ptr<GlobalFrictionData<dims>>& frictionData() const {
- return frictionData_;
- }
-};
#endif
diff --git a/src/data-structures/body.cc b/src/data-structures/body.cc
index b72f396b..39220bd5 100644
--- a/src/data-structures/body.cc
+++ b/src/data-structures/body.cc
@@ -52,36 +52,10 @@ void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assemble() {
assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
- // assemble state energy norm
- typename FrictionBoundaryCondition::BoundaryPatch frictionBoundary(this->leafView(), false); // boundary patch containing all friction boundaries (possibly disconnected)
- for (size_t i=0; i<frictionBoundaryConditions_.size(); i++) {
- const auto& frictionCondition = frictionBoundaryConditions_[i];
- //std::cout << "grid Pointer " << frictionCondition->boundaryPatch()-> << std::endl;
- frictionBoundary.addPatch(*frictionCondition->boundaryPatch());
- }
-
- ScalarMatrix relativeFrictionBoundaryMass;
- assembler_->assembleFrictionalBoundaryMass(frictionBoundary, relativeFrictionBoundaryMass);
- relativeFrictionBoundaryMass /= frictionBoundary.area(); // TODO: weight by individual friction patches?
- stateEnergyNorm_ = std::make_shared<const EnergyNorm<ScalarMatrix, ScalarVector>>(relativeFrictionBoundaryMass);
-
// assemble forces
assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembleFriction(
- const Config::FrictionModel& frictionModel,
- const ScalarVector& weightedNormalStress) {
-
- globalFriction_.resize({frictionBoundaryConditions_.size()});
-
- for (size_t i=0; i<globalFriction_.size(); i++) {
- const auto& frictionBoundaryCondition = frictionBoundaryConditions_[i];
- globalFriction_[i] = assembler_->assembleFrictionNonlinearity(frictionModel, *frictionBoundaryCondition->boundaryPatch(), *frictionBoundaryCondition->frictionData(), weightedNormalStress);
- }
-}
-
// setter and getter
template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
@@ -148,13 +122,6 @@ auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::externalForce() const
return externalForce_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::stateEnergyNorm() const
--> const std::shared_ptr<const StateEnergyNorm>& {
-
- return stateEnergyNorm_;
-}
-
template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition) {
leafBoundaryConditions_.emplace_back(leafBoundaryCondition);
@@ -165,18 +132,10 @@ void Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafBoundaryConditions(
const std::string& tag,
std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const {
- if (tag == "friction") {
- selectedConditions.resize(frictionBoundaryConditions_.size());
- for (size_t i=0; i<frictionBoundaryConditions_.size(); i++) {
- selectedConditions[i] = frictionBoundaryConditions_[i];
- }
-
- } else {
- selectedConditions.resize(0);
- for (size_t i=0; i<leafBoundaryConditions_.size(); i++) {
- if (leafBoundaryConditions_[i]->tag() == tag)
- selectedConditions.emplace_back(leafBoundaryConditions_[i]);
- }
+ selectedConditions.resize(0);
+ for (size_t i=0; i<leafBoundaryConditions_.size(); i++) {
+ if (leafBoundaryConditions_[i]->tag() == tag)
+ selectedConditions.emplace_back(leafBoundaryConditions_[i]);
}
}
@@ -240,26 +199,6 @@ auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions() const
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(std::shared_ptr<FrictionBoundaryCondition> frictionBoundaryCondition) {
- frictionBoundaryConditions_.emplace_back(frictionBoundaryCondition);
-}
-
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::frictionBoundaryConditions() const
--> const std::vector<std::shared_ptr<FrictionBoundaryCondition>>& {
-
- return frictionBoundaryConditions_;
-}
-
-
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::globalFriction()
--> GlobalFrictionContainer& {
- return globalFriction_;
-}
-
-
template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryPatchNodes(
const std::string& tag,
diff --git a/src/data-structures/body.hh b/src/data-structures/body.hh
index 6a254e4c..ed71b586 100644
--- a/src/data-structures/body.hh
+++ b/src/data-structures/body.hh
@@ -24,12 +24,9 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/tectonic/bodydata.hh>
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/globalfrictiondata.hh>
#include "../assemblers.hh"
#include "../boundarycondition.hh"
-#include "globalfrictioncontainer.hh"
#include "enums.hh"
#include "matrices.hh"
@@ -56,8 +53,7 @@ class Body {
using Function = Dune::VirtualFunction<double, double>;
using ExternalForce = std::function<void(double, VectorType &)>;
- using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
- using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
+ using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
using LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
using BoundaryFunctions = std::vector<const Function* >;
@@ -65,9 +61,6 @@ class Body {
using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
using BoundaryPatches = std::vector<const typename LeafBoundaryCondition::BoundaryPatch* >;
- using GlobalFriction = GlobalFriction<Matrix, VectorType>;
- using GlobalFrictionContainer = GlobalFrictionContainer<GlobalFriction, 1>;
-
using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
private:
@@ -89,25 +82,17 @@ class Body {
ExternalForce externalForce_;
VectorType gravityFunctional_;
- std::shared_ptr<const StateEnergyNorm> stateEnergyNorm_;
// boundary conditions
std::vector<std::shared_ptr<LeafBoundaryCondition>> leafBoundaryConditions_;
std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>> levelBoundaryConditions_; // first index: level, second index: bc on lvl
- // friction boundary conditions
- std::vector<std::shared_ptr<FrictionBoundaryCondition>> frictionBoundaryConditions_;
- GlobalFrictionContainer globalFriction_;
-
public:
Body(
const std::shared_ptr<BodyData<dims>>& bodyData,
const std::shared_ptr<GridType>& grid);
void assemble();
- void assembleFriction(
- const Config::FrictionModel& frictionModel,
- const ScalarVector& weightedNormalStress);
// setter and getter
auto data() const -> const std::shared_ptr<BodyData<dims>>&;
@@ -124,8 +109,6 @@ class Body {
auto externalForce() const -> const ExternalForce&;
- auto stateEnergyNorm() const -> const std::shared_ptr<const StateEnergyNorm>&;
-
void addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition);
void leafBoundaryConditions(
const std::string& tag,
@@ -143,11 +126,6 @@ class Body {
auto levelBoundaryConditions(int level) const -> const std::vector<std::shared_ptr<LevelBoundaryCondition>>&;
auto levelBoundaryConditions() const -> const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>&;
- void addBoundaryCondition(std::shared_ptr<FrictionBoundaryCondition> frictionBoundaryCondition);
- auto frictionBoundaryConditions() const -> const std::vector<std::shared_ptr<FrictionBoundaryCondition>>&;
-
- auto globalFriction() -> GlobalFrictionContainer&;
-
void boundaryPatchNodes(
const std::string& tag,
BoundaryPatchNodes& nodes) const;
diff --git a/src/data-structures/levelcontactnetwork.cc b/src/data-structures/levelcontactnetwork.cc
index bf3ca676..9709a867 100644
--- a/src/data-structures/levelcontactnetwork.cc
+++ b/src/data-structures/levelcontactnetwork.cc
@@ -4,8 +4,14 @@
#include <dune/contact/projections/normalprojection.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/tectonic/globalratestatefriction.hh>
+
#include "levelcontactnetwork.hh"
+#include "../assemblers.hh"
#include "../utils/tobool.hh"
#include "../utils/debugutils.hh"
@@ -24,22 +30,50 @@ LevelContactNetwork<GridType, dimension>::LevelContactNetwork(
nBodyAssembler_(nBodies, nCouplings)
{}
+template <class GridType, int dimension>
+LevelContactNetwork<GridType, dimension>::~LevelContactNetwork() {
+ for (size_t i=0; i<nBodies(); i++) {
+ delete frictionBoundaries_[i];
+ delete stateEnergyNorms_[i];
+ }
+}
+
template <class GridType, int dimension>
void LevelContactNetwork<GridType, dimension>::assemble() {
+ using BoundaryPatch = DeformedLeafBoundaryPatch;
+
for (size_t i=0; i<nBodies(); i++) {
bodies_[i]->assemble();
//printBasisDofLocation(bodies_[i]->assembler()->vertexBasis); //TODO remove after debugging
+
+ frictionBoundaries_[i] = new BoundaryPatch(bodies_[i]->leafView(),false);
}
// set up dune-contact nBodyAssembler
nBodyAssembler_.setGrids(deformedGrids_);
for (size_t i=0; i<nCouplings(); i++) {
- nBodyAssembler_.setCoupling(*couplings_[i], i);
+ auto& coupling = couplings_[i];
+ nBodyAssembler_.setCoupling(*coupling, i);
+
+ const auto& contactCoupling = nBodyAssembler_.getContactCouplings()[i]; // dual mortar object holding boundary patches
+ const auto nonmortarIdx = coupling->gridIdx_[0];
+ const auto mortarIdx = coupling->gridIdx_[1];
+
+ frictionBoundaries_[nonmortarIdx]->addPatch(contactCoupling->nonmortarBoundary());
+ frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
}
nBodyAssembler_.assembleTransferOperator();
nBodyAssembler_.assembleObstacle();
+
+ // assemble state energy norm
+ for (size_t i=0; i<nBodies(); i++) {
+ ScalarMatrix relativeFrictionBoundaryMass;
+ bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], relativeFrictionBoundaryMass);
+ relativeFrictionBoundaryMass /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
+ stateEnergyNorms_[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(relativeFrictionBoundaryMass);
+ }
}
template <class GridType, int dimension>
@@ -47,13 +81,39 @@ void LevelContactNetwork<GridType, dimension>::assembleFriction(
const Config::FrictionModel& frictionModel,
const std::vector<ScalarVector>& weightedNormalStress) {
- for (size_t i=0; i<nBodies(); i++) {
- bodies_[i]->assembleFriction(frictionModel, weightedNormalStress[i]);
- }
+ assert(weightedNormalStress.size() == bodies_.size());
+
+ const size_t nBodies = bodies_.size();
+
+ // Lumping of the nonlinearity
+ std::vector<ScalarVector> weights(nBodies);
+ {
+ NeumannBoundaryAssembler<DeformedGridType, typename ScalarVector::block_type>
+ frictionalBoundaryAssembler(std::make_shared<
+ ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+ 1));
+
+ for (size_t i=0; i<nBodies; i++) {
+ bodies_[i]->assembler()->assembleBoundaryFunctional(frictionalBoundaryAssembler, weights[i], *frictionBoundaries_[i]);
+ }
+ }
+
+ switch (frictionModel) {
+ case Config::Truncated:
+ globalFriction_ = std::make_shared<GlobalRateStateFriction<
+ Matrix, Vector, TruncatedRateState, DeformedGridType>>(
+ nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+ case Config::Regularised:
+ globalFriction_ = std::make_shared<GlobalRateStateFriction<
+ Matrix, Vector, RegularisedRateState, DeformedGridType>>(
+ nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+ default:
+ assert(false);
+ }
}
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::shared_ptr<Body>>& bodies) {
+void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::shared_ptr<Body>> bodies) {
assert(nBodies()==bodies.size());
bodies_ = bodies;
@@ -80,7 +140,7 @@ void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::
}
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::setCouplings(const std::vector<std::shared_ptr<CouplingPair>>& couplings) {
+void LevelContactNetwork<GridType, dimension>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
assert(this->nCouplings()==couplings.size());
couplings_ = couplings;
}
@@ -99,13 +159,17 @@ void LevelContactNetwork<GridType, dimension>::constructCoupling(
int nonMortarBodyIdx, int mortarBodyIdx,
const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
- std::shared_ptr<CouplingPair>& coupling) const {
+ const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+ const std::shared_ptr<GlobalFrictionData<dimension>>& frictionData,
+ std::shared_ptr<FrictionCouplingPair>& coupling) const {
- coupling = std::make_shared<CouplingPair>();
+ coupling = std::make_shared<FrictionCouplingPair>();
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<DeformedLeafBoundaryPatch>>();
- std::shared_ptr<typename CouplingPair::BackEndType> backend = nullptr;
+ std::shared_ptr<typename FrictionCouplingPair::BackEndType> backend = nullptr;
coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
+ coupling->setWeakeningPatch(weakeningPatch);
+ coupling->setFrictionData(frictionData);
}
// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
@@ -191,13 +255,20 @@ void LevelContactNetwork<GridType, dimension>::boundaryPatches(
}
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::stateEnergyNorms(StateEnergyNorms& stateEnergyNorms) const {
- stateEnergyNorms.resize(nBodies());
- for (size_t i=0; i<stateEnergyNorms.size(); i++) {
- stateEnergyNorms[i] = body(i)->stateEnergyNorm().get();
+void LevelContactNetwork<GridType, dimension>::frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const {
+ nodes.resize(nBodies());
+ for (size_t i=0; i<nBodies(); i++) {
+ nodes[i] = frictionBoundaries_[i]->getVertices();
}
}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::stateEnergyNorms() const
+-> const StateEnergyNorms& {
+ return stateEnergyNorms_;
+}
+
template <class GridType, int dimension>
auto LevelContactNetwork<GridType, dimension>::level() const
-> int {
@@ -233,6 +304,13 @@ auto LevelContactNetwork<GridType, dimension>::body(int i) const
return bodies_[i];
}
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::coupling(int i) const
+-> const std::shared_ptr<FrictionCouplingPair>& {
+
+ return couplings_[i];
+}
+
template <class GridType, int dimension>
auto LevelContactNetwork<GridType, dimension>::deformedGrids() const
-> const std::vector<const DeformedGridType*>& {
@@ -240,6 +318,13 @@ auto LevelContactNetwork<GridType, dimension>::deformedGrids() const
return deformedGrids_;
}
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::couplings() const
+-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
+
+ return couplings_;
+}
+
template <class GridType, int dimension>
auto LevelContactNetwork<GridType, dimension>::leafView(int bodyID) const
-> const DeformedLeafGridView& {
@@ -290,12 +375,10 @@ auto LevelContactNetwork<GridType, dimension>::matrices() const
}
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::globalFriction(GlobalFrictionContainer& globalFriction) const {
- globalFriction.resize({nBodies()});
+auto LevelContactNetwork<GridType, dimension>::globalFriction() const
+-> GlobalFriction& {
- for (size_t i=0; i<nBodies(); i++) {
- globalFriction[i] = bodies_[i]->globalFriction();
- }
+ return *globalFriction_;
}
template <class GridType, int dimension>
@@ -307,5 +390,4 @@ void LevelContactNetwork<GridType, dimension>::externalForces(ExternalForces& ex
}
}
-
#include "levelcontactnetwork_tmpl.cc"
diff --git a/src/data-structures/levelcontactnetwork.hh b/src/data-structures/levelcontactnetwork.hh
index e7672cd5..4721208b 100644
--- a/src/data-structures/levelcontactnetwork.hh
+++ b/src/data-structures/levelcontactnetwork.hh
@@ -13,11 +13,13 @@
#include <dune/contact/projections/normalprojection.hh>
#include <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
#include "../assemblers.hh"
+#include "../frictioncouplingpair.hh"
#include "body.hh"
#include "enums.hh"
-#include "globalfrictioncontainer.hh"
#include "matrices.hh"
//#include "multi-body-problem-data/myglobalfrictiondata.hh"
@@ -42,7 +44,7 @@ template <class GridType, int dimension> class LevelContactNetwork {
using ScalarVector = typename Assembler::ScalarVector;
using field_type = typename Matrix::field_type;
- using CouplingPair = Dune::Contact::CouplingPair<DeformedGridType, DeformedGridType, field_type>;
+ using FrictionCouplingPair = FrictionCouplingPair<DeformedGridType, LocalVector, field_type>;
using Function = Dune::VirtualFunction<double, double>;
@@ -51,7 +53,7 @@ template <class GridType, int dimension> class LevelContactNetwork {
using BoundaryPatchNodes = std::vector<typename Body::BoundaryPatchNodes>;
using BoundaryPatches = std::vector<typename Body::BoundaryPatches>;
- using GlobalFrictionContainer = GlobalFrictionContainer<typename Body::GlobalFriction, 2>;
+ using GlobalFriction = GlobalFriction<Matrix, Vector>;
using StateEnergyNorms = std::vector<const typename Body::StateEnergyNorm *>;
using ExternalForces = std::vector<const typename Body::ExternalForce *>;
@@ -60,14 +62,15 @@ template <class GridType, int dimension> class LevelContactNetwork {
public:
LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
+ ~LevelContactNetwork();
void assemble();
void assembleFriction(
const Config::FrictionModel& frictionModel,
const std::vector<ScalarVector>& weightedNormalStress);
- void setBodies(const std::vector<std::shared_ptr<Body>>& bodies);
- void setCouplings(const std::vector<std::shared_ptr<CouplingPair>>& couplings);
+ void setBodies(const std::vector<std::shared_ptr<Body>> bodies);
+ void setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings);
void constructBody(
const std::shared_ptr<BodyData<dimension>>& bodyData,
@@ -78,7 +81,9 @@ template <class GridType, int dimension> class LevelContactNetwork {
int mortarBodyIdx,
const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
- std::shared_ptr<CouplingPair>& coupling) const;
+ const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+ const std::shared_ptr<GlobalFrictionData<dimension>>& frictionData,
+ std::shared_ptr<FrictionCouplingPair>& coupling) const;
// getter
void totalNodes(
@@ -96,9 +101,10 @@ template <class GridType, int dimension> class LevelContactNetwork {
void boundaryFunctions(
const std::string& tag,
BoundaryFunctions& functions) const;
+ void frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const;
- void stateEnergyNorms(StateEnergyNorms& stateEnergyNorms) const;
- void globalFriction(GlobalFrictionContainer& globalFriction) const;
+ auto stateEnergyNorms() const -> const StateEnergyNorms&;
+ auto globalFriction() const -> GlobalFriction&;
void externalForces(ExternalForces& externalForces) const;
auto level() const -> int;
@@ -108,6 +114,9 @@ template <class GridType, int dimension> class LevelContactNetwork {
auto body(int i) const -> const std::shared_ptr<Body>&;
+ auto coupling(int i) const -> const std::shared_ptr<FrictionCouplingPair>&;
+ auto couplings() const -> const std::vector<std::shared_ptr<FrictionCouplingPair>>&;
+
auto deformedGrids() const -> const std::vector<const DeformedGridType*>& ;
auto leafVertexCounts() const -> const std::vector<size_t>&;
@@ -122,21 +131,24 @@ template <class GridType, int dimension> class LevelContactNetwork {
auto matrices() const -> const Matrices<Matrix,2>&;
-
-
private:
const int level_;
std::vector<std::shared_ptr<Body>> bodies_;
- std::vector<std::shared_ptr<CouplingPair>> couplings_;
+ std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
std::vector<const DeformedGridType*> deformedGrids_;
std::vector<std::unique_ptr<const DeformedLeafGridView>> leafViews_;
std::vector<std::vector<std::unique_ptr<const DeformedLevelGridView>>> levelViews_;
std::vector<size_t> leafVertexCounts_;
+ std::vector<DeformedLeafBoundaryPatch*> frictionBoundaries_; // union of all boundary patches per body
+ StateEnergyNorms stateEnergyNorms_;
+
NBodyAssembler nBodyAssembler_;
Matrices<Matrix,2> matrices_;
+
+ std::shared_ptr<GlobalFriction> globalFriction_;
};
#endif
diff --git a/src/data-structures/levelcontactnetwork_tmpl.cc b/src/data-structures/levelcontactnetwork_tmpl.cc
index 1a4be09f..1bd54512 100644
--- a/src/data-structures/levelcontactnetwork_tmpl.cc
+++ b/src/data-structures/levelcontactnetwork_tmpl.cc
@@ -5,4 +5,6 @@
#include "../explicitgrid.hh"
#include "levelcontactnetwork.hh"
+using MyLevelContactNetwork = LevelContactNetwork<Grid, MY_DIM>;
+
template class LevelContactNetwork<Grid, MY_DIM>;
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 181e5c86..ec0311fe 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -151,7 +151,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
print(upper, "upper");
// set up functional and TNMMG solver
- using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity, Vector&, Vector&, typename Matrix::field_type>;
+ using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
diff --git a/src/factories/levelcontactnetworkfactory.hh b/src/factories/levelcontactnetworkfactory.hh
index fa026527..833036c6 100644
--- a/src/factories/levelcontactnetworkfactory.hh
+++ b/src/factories/levelcontactnetworkfactory.hh
@@ -12,14 +12,14 @@ class LevelContactNetworkFactory {
protected:
using Body = typename LevelContactNetwork::Body;
- using CouplingPair = typename LevelContactNetwork::CouplingPair;
+ using FrictionCouplingPair = typename LevelContactNetwork::FrictionCouplingPair;
const Dune::ParameterTree& parset_;
const size_t bodyCount_;
const size_t couplingCount_;
std::vector<std::shared_ptr<Body>> bodies_;
- std::vector<std::shared_ptr<CouplingPair>> couplings_;
+ std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
LevelContactNetwork levelContactNetwork_;
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index 11a3f98c..0e161498 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -114,22 +114,23 @@ void StackedBlocksFactory<GridType, dims>::setCouplings() {
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
- std::shared_ptr<typename Base::CouplingPair::BackEndType> backend = nullptr;
+ std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
for (size_t i=0; i<this->couplings_.size(); i++) {
auto& coupling = this->couplings_[i];
- coupling = std::make_shared<typename Base::CouplingPair>();
+ coupling = std::make_shared<typename Base::FrictionCouplingPair>();
nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
- coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::CouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+ coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+ coupling->setWeakeningPatch(upperWeakPatches_[i]);
+ coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], CuboidGeometry::lengthScale));
}
}
template <class GridType, int dims>
void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
- using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
@@ -143,22 +144,6 @@ void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
- // friction boundary
- if (i<this->bodyCount_-1 && upperWeakPatches_[i]->vertices.size()>0) {
- std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
- frictionBoundary->setBoundaryPatch(leafFaces_[i]->upper);
- frictionBoundary->setWeakeningPatch(upperWeakPatches_[i]);
- frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], lengthScale));
- body->addBoundaryCondition(frictionBoundary);
- }
- if (i>0 && lowerWeakPatches_[i]->vertices.size()>0) {
- std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
- frictionBoundary->setBoundaryPatch(leafFaces_[i]->lower);
- frictionBoundary->setWeakeningPatch(lowerWeakPatches_[i]);
- frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *lowerWeakPatches_[i], lengthScale));
- body->addBoundaryCondition(frictionBoundary);
- }
-
// Neumann boundary
std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
body->addBoundaryCondition(neumannBoundary);
diff --git a/src/frictioncouplingpair.hh b/src/frictioncouplingpair.hh
new file mode 100644
index 00000000..7eb36871
--- /dev/null
+++ b/src/frictioncouplingpair.hh
@@ -0,0 +1,37 @@
+#ifndef SRC_FRICTIONCOUPLINGPAIR_HH
+#define SRC_FRICTIONCOUPLINGPAIR_HH
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+template <class GridType, class LocalVectorType, class field_type = double>
+class FrictionCouplingPair : public Dune::Contact::CouplingPair<GridType, GridType, field_type>{
+private:
+ static const int dim = GridType::dimensionworld;
+
+ using Base = Dune::Contact::CouplingPair<GridType,GridType,field_type>;
+ using LocalVector = LocalVectorType;
+
+ // friction data
+ std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
+ std::shared_ptr<GlobalFrictionData<dim>> frictionData_;
+
+public:
+ void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
+ weakeningPatch_ = weakeningPatch;
+ }
+
+ void setFrictionData(std::shared_ptr<GlobalFrictionData<dim>> frictionData) {
+ frictionData_ = frictionData;
+ }
+
+ const auto& weakeningPatch() const {
+ return *weakeningPatch_;
+ }
+
+ const auto& frictionData() const {
+ return *frictionData_;
+ }
+};
+#endif
diff --git a/src/io/hdf5/frictionalboundary-writer.cc b/src/io/hdf5/frictionalboundary-writer.cc
index 968b404f..19641aa5 100644
--- a/src/io/hdf5/frictionalboundary-writer.cc
+++ b/src/io/hdf5/frictionalboundary-writer.cc
@@ -52,7 +52,7 @@ void FrictionalBoundaryWriter<ProgramState, GridView>::write(
addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
frictionalBoundaryStates);
- friction.updateAlpha(programState.alpha[id_]);
+ friction.updateAlpha(programState.alpha);
auto const c = friction.coefficientOfFriction(programState.v[id_]);
auto const frictionalBoundaryCoefficient =
restrictToSurface(c, frictionalBoundary_);
diff --git a/src/multi-body-problem-data/weakpatch.hh b/src/multi-body-problem-data/weakpatch.hh
index a1fde1cb..047dc482 100644
--- a/src/multi-body-problem-data/weakpatch.hh
+++ b/src/multi-body-problem-data/weakpatch.hh
@@ -10,7 +10,7 @@
template <class LocalVector>
void getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry, ConvexPolyhedron<LocalVector>& lowerWeakPatch, ConvexPolyhedron<LocalVector>& upperWeakPatch) {
-#if MY_DIM == 3
+#if MY_DIM == 3 // TODO: Does not work yet
if (cuboidGeometry.V != cuboidGeometry.W) {
lowerWeakPatch.vertices.resize(4);
lowerWeakPatch.vertices[0] = lowerWeakPatch.vertices[2] = cuboidGeometry.V;
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 4d6e7173..709adff5 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -192,8 +192,7 @@ int main(int argc, char *argv[]) {
// ------------------------
levelContactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
- typename LevelContactNetwork::GlobalFrictionContainer globalFriction;
- levelContactNetwork.globalFriction(globalFriction);
+ auto& globalFriction = levelContactNetwork.globalFriction();
globalFriction.updateAlpha(programState.alpha);
using MyVertexBasis = typename Assembler::VertexBasis;
@@ -285,17 +284,14 @@ int main(int argc, char *argv[]) {
using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
StateUpdater<ScalarVector, Vector>>;
- std::vector<double> L(bodyCount, parset.get<double>("boundary.friction.L"));
- std::vector<double> V0(bodyCount, parset.get<double>("boundary.friction.V0"));
-
BoundaryFunctions velocityDirichletFunctions;
levelContactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
BoundaryNodes dirichletNodes;
levelContactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- typename LevelContactNetwork::BoundaryPatchNodes frictionNodes;
- levelContactNetwork.boundaryPatchNodes("friction", frictionNodes);
+ std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+ levelContactNetwork.frictionNodes(frictionNodes);
Updaters current(
initRateUpdater(
@@ -306,18 +302,17 @@ int main(int argc, char *argv[]) {
programState.u,
programState.v,
programState.a),
- initStateUpdater<ScalarVector, Vector, typename LevelContactNetwork::BoundaryPatchNodes>(
+ initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
programState.alpha,
- frictionNodes,
- L,
- V0));
+ nBodyAssembler.getContactCouplings(),
+ levelContactNetwork.couplings())
+ );
auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
- typename LevelContactNetwork::StateEnergyNorms stateEnergyNorms;
- levelContactNetwork.stateEnergyNorms(stateEnergyNorms);
+ const auto& stateEnergyNorms = levelContactNetwork.stateEnergyNorms();
auto const mustRefine = [&](Updaters &coarseUpdater,
Updaters &fineUpdater) {
@@ -342,8 +337,8 @@ int main(int argc, char *argv[]) {
typename LevelContactNetwork::ExternalForces externalForces;
levelContactNetwork.externalForces(externalForces);
- AdaptiveTimeStepper<NonlinearFactory, decltype(nBodyAssembler), Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
- adaptiveTimeStepper(parset, nBodyAssembler, globalFriction, current,
+ AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(nBodyAssembler)>, Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
+ adaptiveTimeStepper(parset, nBodyAssembler, globalFriction, frictionNodes, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 1dbaec37..5ca84926 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -39,10 +39,13 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIterator(
Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
- GlobalFrictionContainer& globalFriction, const std::vector<const ErrorNorm* >& errorNorms)
+ GlobalFriction& globalFriction,
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+ const std::vector<const ErrorNorm* >& errorNorms)
: parset_(parset),
nBodyAssembler_(nBodyAssembler),
globalFriction_(globalFriction),
+ bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
lambda_(parset.get<double>("v.fpi.lambda")),
@@ -57,8 +60,10 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
- std::vector<const Matrix*> matrices_ptr(velocityMatrices.size());
- for (size_t i=0; i<matrices_ptr.size(); i++) {
+ const auto nBodies = nBodyAssembler_.nGrids();
+
+ std::vector<const Matrix*> matrices_ptr(nBodies);
+ for (size_t i=0; i<nBodies; i++) {
matrices_ptr[i] = &velocityMatrices[i];
}
@@ -88,8 +93,8 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
}
// set up functional and TNMMG solver
- Functional J(bilinearForm, totalRhs, ..., lower, upper);
- Factory solverFactory(parset.sub("solver.tnnmg"), J);
+ Functional J(bilinearForm, totalRhs, globalFriction_, lower, upper);
+ Factory solverFactory(parset_.sub("solver.tnnmg"), J);
auto step = solverFactory.step();
EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
@@ -99,50 +104,19 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
size_t fixedPointIteration;
size_t multigridIterations = 0;
- std::vector<ScalarVector> alpha;
+ std::vector<ScalarVector> alpha(nBodies);
updaters.state_->extractAlpha(alpha);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
- // compute relative velocities
- std::vector<Vector> v_rel;
- relativeVelocities(velocityIterates, v_rel);
-
// contribution from nonlinearity
globalFriction_.updateAlpha(alpha);
-
- std::vector<Matrix> matrices(velocityMatrices.size());
- std::vector<Vector> rhs(velocityRHSs.size());
- for (size_t i=0; i<globalFriction_.size(); i++) {
- matrices[i] = velocityMatrices[i];
- rhs[i] = velocityRHSs[i];
-
- auto& globalFriction = globalFriction_[i];
- for (size_t j=0; j<globalFriction.size(); j++) {
- globalFriction[j]->addHessian(v_rel[i], matrices[i]);
- globalFriction[j]->addGradient(v_rel[i], rhs[i]);
- }
-
- matrices_ptr[i] = &matrices[i];
- }
-
- // assemble full global contact problem
- Matrix bilinearForm;
- nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
-
- Vector totalRhs;
- nBodyAssembler_.assembleRightHandSide(rhs, totalRhs);
-
Vector totalVelocityIterate;
nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
- print(bilinearForm, "matrix: ");
- print(totalRhs, "totalRhs: ");
- print(totalVelocityIterate, "iterate: ");
-
// solve a velocity problem
- step->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
+ step->setProblem(totalVelocityIterate);
velocityProblemSolver.preprocess();
velocityProblemSolver.solve();
@@ -158,16 +132,17 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
Dune::MatrixVector::addProduct(v_m[i], lambda_, velocityIterates[i]);
}
- // compute relative velocities on contact boundaries
- relativeVelocities(v_m, v_rel);
+ // extract relative velocities in mortar basis
+ std::vector<Vector> v_rel;
+ relativeVelocities(totalVelocityIterate, v_rel);
// solve a state problem
updaters.state_->solve(v_rel);
- std::vector<ScalarVector> newAlpha;
+ std::vector<ScalarVector> newAlpha(nBodies);
updaters.state_->extractAlpha(newAlpha);
bool breakCriterion = true;
- for (size_t i=0; i<alpha.size(); i++) {
+ for (size_t i=0; i<nBodies; i++) {
if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
breakCriterion = false;
break;
@@ -180,6 +155,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
}
alpha = newAlpha;
}
+
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
@@ -200,189 +176,44 @@ std::ostream &operator<<(std::ostream &stream,
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const {
- // init result
- v_rel.resize(v.size());
- for (size_t i=0; i<v_rel.size(); i++) {
- v_rel[i].resize(v[i].size());
- v_rel[i] = 0;
- }
-
- // needs assemblers to obtain basis
- /*
- std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
-
- using field_type = typename Factory::Matrix::field_type;
-
-
- // adaptation of DualMortarCoupling::setup()
-
- const size_t dim = DeformedGrid::dimension;
- typedef typename DeformedGrid::LeafGridView GridView;
-
- //cache of local bases
- typedef Dune::PQkLocalFiniteElementCache<typename DeformedGrid::ctype, field_type, dim,1> FiniteElementCache1;
- FiniteElementCache1 cache1;
-
- // cache for the dual functions on the boundary
- using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
- std::unique_ptr<DualCache> dualCache;
- dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
-
- // define FE grid functions
- std::vector<BasisGridFunction<VertexBasis, Vector>* > gridFunctions(v_m.size());
- for (size_t i=0; i<gridFunctions.size(); i++) {
- gridFunctions[i] = new BasisGridFunction<MyAssembler::VertexBasis, Vector>(assemblers[i]->vertexBasis, v_m[i]);
- }
- */
-
- /*
- for (size_t i=0; i<nBodyAssembler_.nCouplings(); i++) {
- const auto& coupling = nBodyAssembler_.getCoupling(i);
- auto glue = coupling.backend();
-
- const std::array<int, 2> gridIdx = coupling.gridIdx_;
- const int nonmortarGridIdx = ;
- const int mortarGridIdx = ;
-
- // loop over all intersections
- 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);
+void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeVelocities(
+ const Vector& v,
+ std::vector<Vector>& v_rel) const {
- const auto& mortarFiniteElement = cache1.get(mortarEType);
- dualCache->bind(inside, rIs.indexInInside());
+ const size_t nBodies = nBodyAssembler_.nGrids();
+ // const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
- 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();
+ v_rel.resize(nBodies);
- 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];
+ /* for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+ boundaryNodes =
- weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
- * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+ const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
- // loop over all mortar shape functions
- for (int k=0; k<noOfMortarVec; k++) {
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
- int colIdx = indexSet1.subIndex(outside, k, dim);
- if (!mortarBoundary_.containsVertex(colIdx))
- continue;
+ for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
+ const auto i = vertexMapper.index(*it);
- // Integrate over the product of two shape functions
- field_type mortarEntry = integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+ for (size_t j=0; j<couplingIndices.size(); j++) {
+ const auto couplingIdx = couplingIndices[j];
- Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
-
- }
+ if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
+ continue;
+ localToGlobal_.emplace_back(i);
+ restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
+ couplings[i]->frictionData()(geoToPoint(it->geometry())));
+ break;
}
+ globalIdx++;
}
-
-
-
-
-
-
- }
-
-
- // 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);
- }
- }
- }
- */
-
+ maxIndex_[bodyIdx] = globalIdx;
+ }*/
}
+
#include "fixedpointiterator_tmpl.cc"
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index 02ad3089..078b1d20 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -4,13 +4,14 @@
#include <memory>
#include <dune/common/parametertree.hh>
+#include <dune/common/bitsetvector.hh>
#include <dune/solvers/norms/norm.hh>
#include <dune/solvers/solvers/solver.hh>
-#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
-#include "../data-structures/globalfrictioncontainer.hh"
+#include <dune/contact/assemblers/nbodyassembler.hh>
struct FixedPointIterationCounter {
size_t iterations = 0;
@@ -36,15 +37,17 @@ class FixedPointIterator {
// using DeformedGrid = typename Factory::DeformedGrid;
public:
- using GlobalFrictionContainer = GlobalFrictionContainer<Nonlinearity, 2>;
+ using GlobalFriction = Nonlinearity;
+ using BitVector = Dune::BitSetVector<1>;
private:
- void relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const;
+ void relativeVelocities(const Vector& v, std::vector<Vector>& v_rel) const;
public:
FixedPointIterator(const Dune::ParameterTree& parset,
NBodyAssembler& nBodyAssembler,
- GlobalFrictionContainer& globalFriction,
+ GlobalFriction& globalFriction,
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
const std::vector<const ErrorNorm* >& errorNorms);
FixedPointIterationCounter run(Updaters updaters,
@@ -53,9 +56,11 @@ class FixedPointIterator {
std::vector<Vector>& velocityIterates);
private:
- Dune::ParameterTree const &parset_;
+ const Dune::ParameterTree& parset_;
NBodyAssembler& nBodyAssembler_;
- GlobalFrictionContainer& globalFriction_;
+ GlobalFriction& globalFriction_;
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
size_t fixedPointMaxIterations_;
double fixedPointTolerance_;
diff --git a/src/spatial-solving/fixedpointiterator_tmpl.cc b/src/spatial-solving/fixedpointiterator_tmpl.cc
index 20322c8c..a91a7c6a 100644
--- a/src/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/src/spatial-solving/fixedpointiterator_tmpl.cc
@@ -5,45 +5,24 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
-#include <dune/common/function.hh>
-
#include <dune/solvers/norms/energynorm.hh>
+#include "../spatial-solving/solverfactory_tmpl.cc"
+#include "../data-structures/levelcontactnetwork_tmpl.cc"
-#include <dune/contact/common/deformedcontinuacomplex.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-
-#include "../data-structures/levelcontactnetwork.hh"
#include "../time-stepping/rate/rateupdater.hh"
#include "../time-stepping/state/stateupdater.hh"
#include "../time-stepping/updaters.hh"
-#include "solverfactory.hh"
-
-#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-//#include "../../dune/tectonic/globalfriction.hh"
-//#include "functional.hh"
-//#include "localbisectionsolver.hh"
-
-using MyLevelContactNetwork = LevelContactNetwork<Grid, MY_DIM>;
+using NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
-using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
-using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
-using NonlinearLocalSolver = LocalBisectionSolver;
-using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<MyFunctional, NonlinearLocalSolver>;
-using Factory = SolverFactory<MyFunctional, BitVector>;
-
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-template class FixedPointIterator<Factory, typename MyLevelContactNetwork::NBodyAssembler, MyUpdaters, ErrorNorm>;
+
+template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 906649b6..c9a3d3b1 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -14,7 +14,7 @@
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
#include <dune/tnnmg/projections/obstacledefectprojection.hh>
-#include "tnnmg/functional.hh"
+//#include "tnnmg/functional.hh"
#include "tnnmg/linearization.hh"
#include "tnnmg/linesearchsolver.hh"
#include "tnnmg/localbisectionsolver.hh"
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 60fc8802..1355a1fe 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -5,10 +5,15 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
-//#include "../../dune/tectonic/globalratestatefriction.hh"
+#include "../../dune/tectonic/globalfriction.hh"
#include "tnnmg/functional.hh"
#include "tnnmg/zerononlinearity.hh"
+#include "solverfactory.hh"
-using FrictionlessFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity, Vector&, Vector&, double>;
-template class SolverFactory<FrictionlessFunctional, BitVector>;
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MyZeroFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+template class SolverFactory<MyFunctional, BitVector>;
+template class SolverFactory<MyZeroFunctional, BitVector>;
diff --git a/src/tests/contactmerge.cc b/src/tests/contactmerge.cc
new file mode 100644
index 00000000..5bf8d918
--- /dev/null
+++ b/src/tests/contactmerge.cc
@@ -0,0 +1,40 @@
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid-glue/adapter/gridgluevtkwriter.hh>
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+#include <dune/grid-glue/gridglue.hh>
+#include <dune/grid-glue/merging/contactmerge.hh>
+
+const unsigned dim = 3;
+using Coordinates = Dune::EquidistantOffsetCoordinates<double, dim>;
+using Grid = Dune::YaspGrid<dim, Coordinates>;
+using Element = Grid::Codim<0>::Entity;
+using Extractor = Dune::GridGlue::Codim1Extractor<Grid::LeafGridView>;
+using GridGlue = Dune::GridGlue::GridGlue<Extractor, Extractor>;
+using ContactMerge = Dune::GridGlue::ContactMerge<dim, Grid::ctype>;
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ Grid grid0{{0., 0., 0.}, {1., 1., 1.}, {10, 10, 10}};
+ Grid grid1{{.12, 0.23, 1.05}, {1.12, 1.23, 2.05}, {10, 10, 10}};
+
+ auto truePredicate = [](const Element&, unsigned int) { return true; };
+
+ auto extractor0 = std::make_shared<Extractor>(grid0.leafGridView(), truePredicate);
+ auto extractor1 = std::make_shared<Extractor>(grid1.leafGridView(), truePredicate);
+
+ auto merger = std::make_shared<ContactMerge>();
+
+ GridGlue glue(extractor0, extractor1, merger);
+ glue.build();
+
+ Dune::GridGlue::GridGlueVtkWriter::write(glue, "contactmerge");
+
+ return 0;
+}
diff --git a/src/tests/couplingtest.hh b/src/tests/couplingtest.hh
new file mode 100644
index 00000000..8346e5dd
--- /dev/null
+++ b/src/tests/couplingtest.hh
@@ -0,0 +1,192 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GRIDGLUE_TEST_COUPLINGTEST_HH
+#define DUNE_GRIDGLUE_TEST_COUPLINGTEST_HH
+
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/version.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/grid-glue/gridglue.hh>
+
+template <class IntersectionIt>
+bool testIntersection(const IntersectionIt & rIIt)
+{
+ bool success = true;
+
+ typedef typename IntersectionIt::value_type Intersection;
+ // Dimension of the intersection
+ const int dim = Intersection::mydim;
+
+ // Dimension of world coordinates
+ const int coorddim = Intersection::coorddim;
+
+ // Create a set of test points
+ const Dune::QuadratureRule<double, dim>& quad = Dune::QuadratureRules<double, dim>::rule(rIIt->type(), 3);
+
+ for (unsigned int l=0; l<quad.size(); l++) {
+ const auto inside = rIIt->inside();
+ const auto outside = rIIt->outside();
+
+ Dune::FieldVector<double, dim> quadPos = quad[l].position();
+
+ Dune::FieldVector<double, Intersection::InsideGridView::dimensionworld> localGrid0Pos =
+ inside.geometry().global(rIIt->geometryInInside().global(quadPos));
+
+ // currently the intersection maps to the GV::dimworld, this will hopefully change soon
+ Dune::FieldVector<double, Intersection::InsideGridView::dimensionworld> globalGrid0Pos =
+ rIIt->geometry().global(quadPos);
+
+ Dune::FieldVector<double, Intersection::OutsideGridView::dimensionworld> localGrid1Pos =
+ outside.geometry().global(rIIt->geometryInOutside().global(quadPos));
+
+ // currently the intersection maps to the GV::dimworld, this will hopefully change soon
+ Dune::FieldVector<double, Intersection::OutsideGridView::dimensionworld> globalGrid1Pos =
+ rIIt->geometryOutside().global(quadPos);
+
+ // Test whether local grid0 position is consistent with global grid0 position
+ if ( (localGrid0Pos-globalGrid0Pos).two_norm() >= 1e-6 )
+ {
+ std::cout << __FILE__ << ":" << __LINE__ << ": error: assert( (localGrid0Pos-globalGrid0Pos).two_norm() < 1e-6 ) failed\n";
+ std::cerr << "localGrid0Pos = " << localGrid0Pos << "\n";
+ std::cerr << "globalGrid0Pos = " << globalGrid0Pos << "\n";
+ success = false;
+ }
+
+ // Test whether local grid1 position is consistent with global grid1 position
+ if ( (localGrid1Pos-globalGrid1Pos).two_norm() >= 1e-6 )
+ {
+ std::cout << __FILE__ << ":" << __LINE__ << ": error: assert( (localGrid1Pos-globalGrid1Pos).two_norm() < 1e-6 ) failed\n";
+ std::cerr << "localGrid1Pos = " << localGrid1Pos << "\n";
+ std::cerr << "globalGrid1Pos = " << globalGrid1Pos << "\n";
+ success = false;
+ }
+
+ // Here we assume that the two interfaces match geometrically:
+ if ( (globalGrid0Pos-globalGrid1Pos).two_norm() >= 1e-4 )
+ {
+ std::cout << __FILE__ << ":" << __LINE__ << ": error: assert( (globalGrid0Pos-globalGrid1Pos).two_norm() < 1e-4 ) failed\n";
+ std::cerr << "localGrid0Pos = " << localGrid0Pos << "\n";
+ std::cerr << "globalGrid0Pos = " << globalGrid0Pos << "\n";
+ std::cerr << "localGrid1Pos = " << localGrid1Pos << "\n";
+ std::cerr << "globalGrid1Pos = " << globalGrid1Pos << "\n";
+ success = false;
+ }
+
+ // Test the normal vector methods. At least test whether they don't crash
+ if (coorddim - dim == 1) // only test for codim 1
+ {
+ rIIt->outerNormal(quadPos);
+ rIIt->unitOuterNormal(quadPos);
+ rIIt->integrationOuterNormal(quadPos);
+ rIIt->centerUnitOuterNormal();
+ }
+ }
+
+ return success;
+}
+
+
+template <class GlueType>
+void testCoupling(const GlueType& glue)
+{
+ bool success = true;
+#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6)
+ using View0Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<0> >;
+ using View1Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<1> >;
+ View0Mapper view0mapper(glue.template gridView<0>(), Dune::mcmgElementLayout());
+ View1Mapper view1mapper(glue.template gridView<1>(), Dune::mcmgElementLayout());
+#else
+ using View0Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<0>, Dune::MCMGElementLayout >;
+ using View1Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<1>, Dune::MCMGElementLayout >;
+ View0Mapper view0mapper(glue.template gridView<0>());
+ View1Mapper view1mapper(glue.template gridView<1>());
+#endif
+
+ std::vector<unsigned int> countInside0(view0mapper.size());
+ std::vector<unsigned int> countOutside1(view1mapper.size());
+ std::vector<unsigned int> countInside1(view1mapper.size(), 0);
+ std::vector<unsigned int> countOutside0(view0mapper.size(), 0);
+
+ // ///////////////////////////////////////
+ // IndexSet
+ // ///////////////////////////////////////
+
+ {
+ size_t count = 0;
+ for (auto rIIt = glue.template ibegin<0>(); rIIt != glue.template iend<0>(); ++rIIt) count ++;
+ typename GlueType::IndexSet is = glue.indexSet();
+ if(is.size() != glue.size())
+ DUNE_THROW(Dune::Exception,
+ "Inconsistent size information: indexSet.size() " << is.size() << " != GridGlue.size() " << glue.size());
+ if(is.size() != count)
+ DUNE_THROW(Dune::Exception,
+ "Inconsistent size information: indexSet.size() " << is.size() << " != iterator count " << count);
+ std::vector<bool> visited(count, false);
+ for (auto rIIt = glue.template ibegin<0>(); rIIt != glue.template iend<0>(); ++rIIt) {
+ size_t idx = is.index(*rIIt);
+ if(idx >= count)
+ DUNE_THROW(Dune::Exception,
+ "Inconsistent IndexSet: index " << idx << " out of range, size is " << count);
+ if(visited[idx] != false)
+ DUNE_THROW(Dune::Exception,
+ "Inconsistent IndexSet: visited index " << idx << " twice");
+ visited[idx] = true;
+ }
+ // make sure that we have a consecutive zero starting index set
+ for (size_t i = 0; i<count; i++)
+ {
+ if (visited[i] != true)
+ DUNE_THROW(Dune::Exception,
+ "Non-consective IndexSet: " << i << " missing.");
+ }
+ }
+
+
+ // ///////////////////////////////////////
+ // MergedGrid centric Grid0->Grid1
+ // ///////////////////////////////////////
+
+ {
+ for (auto rIIt = glue.template ibegin<0>(); rIIt != glue.template iend<0>(); ++rIIt)
+ {
+ if (rIIt->self() && rIIt->neighbor())
+ {
+ const auto index0 = view0mapper.index(rIIt->inside());
+ const auto index1 = view1mapper.index(rIIt->outside());
+
+ countInside0[index0]++;
+ countOutside1[index1]++;
+ success = success && testIntersection(rIIt);
+ }
+ }
+ }
+
+ // ///////////////////////////////////////
+ // MergedGrid centric Grid1->Grid0
+ // ///////////////////////////////////////
+
+ {
+ for (auto rIIt = glue.template ibegin<1>(); rIIt != glue.template iend<1>(); ++rIIt)
+ {
+ if (rIIt->self() && rIIt->neighbor())
+ {
+ const auto index1 = view1mapper.index(rIIt->inside());
+ const auto index0 = view0mapper.index(rIIt->outside());
+
+ countInside1[index1]++;
+ countOutside0[index0]++;
+ success = success && testIntersection(rIIt);
+ }
+ }
+ }
+
+ if (! success)
+ DUNE_THROW(Dune::Exception, "Test failed, see above for details.");
+}
+
+#endif // DUNE_GRIDGLUE_TEST_COUPLINGTEST_HH
diff --git a/src/tests/gridgluefrictiontest.cc b/src/tests/gridgluefrictiontest.cc
index 4e5ee68d..c33a8493 100644
--- a/src/tests/gridgluefrictiontest.cc
+++ b/src/tests/gridgluefrictiontest.cc
@@ -12,11 +12,15 @@
#include <dune/common/stdstreams.hh>
#include <dune/common/fvector.hh>
#include <dune/common/function.hh>
+#include <dune/common/bitsetvector.hh>
//#include <dune/fufem/assemblers/assembler.hh>
#include <dune/fufem/functiontools/basisinterpolator.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
@@ -29,10 +33,11 @@
#include <dune/contact/common/couplingpair.hh>
#include <dune/contact/assemblers/dualmortarcoupling.hh>
#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
#include "../utils/debugutils.hh"
-#include <dune/tectonic/transformedglobalratestatefriction.hh>
+//#include <dune/tectonic/transformedglobalratestatefriction.hh>
#include "common.hh"
@@ -70,9 +75,9 @@ class BoundaryFunction0 : public Dune::VirtualFunction<VectorType, c_type> {
public:
void evaluate(const VectorType& x, c_type& res) const override {
if (isClose(x[1], 1.0)) {
- res = std::sin(x[0]);
+ res = 1.0 ;// +std::sin(x[0]);
} else {
- res = 0.0;
+ res = 1.0;
}
}
};
@@ -83,13 +88,129 @@ class BoundaryFunction1 : public Dune::VirtualFunction<VectorType, c_type> {
public:
void evaluate(const VectorType& x, c_type& res) const override{
if (isClose(x[1], 1.0)) {
- res = std::cos(x[0]);
+ res = 2.0 ;//+ std::cos(x[0]);
} else {
- res = 0.0;
+ res = 2.0;
}
}
};
+/*
+template <class GridView>
+void dualWeights(const BoundaryPatch<GridView>& boundaryPatch, BlockVector& weights, bool initializeVector = true) {
+ typedef typename BoundaryPatch<GridView>::iterator BoundaryIterator;
+
+ /*typedef typename LocalBoundaryFunctionalAssemblerType::LocalVector LocalVector;
+ typedef typename TrialBasis::LinearCombination LinearCombination;
+
+ int rows = tBasis_.size();
+
+
+ const auto inside = it->inside();
+
+ // get shape functions
+ DualP1LocalFiniteElement<class D, class R, dim>
+ class
+ const typename TrialBasis::LocalFiniteElement& lFE = tBasis_.getLocalFiniteElement(inside);
+
+ LocalVector localB(tFE.localBasis().size());
+
+ }*/
+/*
+ // cache for the dual functions on the boundary
+ using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
+ std::unique_ptr<DualCache> dualCache;
+ dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
+
+ if (initializeVector) {
+ weights.resize(rows);
+ std::fill(weights.begin(), weights.end(), 0.0);
+ }
+
+ // loop over all boundary intersections
+ BoundaryIterator it = boundaryPatch.begin();
+ BoundaryIterator end = boundaryPatch.end();
+
+ for (; it != end; ++it) {
+ const auto& inside = it->inside();
+
+ if (!boundaryPatch.contains(inside, it->indexInInside()))
+ continue;
+
+ // types of the elements supporting the boundary segments in question
+ Dune::GeometryType nonmortarEType = inside.type();
+
+ const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+
+ int noOfMortarVec = targetRefElement.size(dim);
+
+ Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),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());
+ const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(it->type(), quadOrder);
+
+ dualCache->bind(inside, it->indexInInside());
+
+
+
+ const auto& rGeomInInside = it->geometryInInside();
+
+ int nNonmortarFaceNodes = domainRefElement.size(it->indexInInside(),1,dim);
+ std::vector<int> nonmortarFaceNodes;
+ for (int i=0; i<nNonmortarFaceNodes; i++) {
+ int faceIdxi = domainRefElement.subEntity(it->indexInInside(), 1, i, dim);
+ nonmortarFaceNodes.push_back(faceIdxi);
+ }
+
+ // store values of shape functions
+ std::vector<Dune::FieldVector<field_type,1>> values(tFE.localBasis().size());
+
+ // get geometry and store it
+ const auto& rGeom = it->geometry();
+
+ localVector = 0.0;
+
+ // get quadrature rule
+ QuadratureRuleKey tFEquad(tFE);
+ QuadratureRuleKey quadKey = tFEquad.product(functionQuadKey_);
+ const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
+
+ // loop over quadrature points
+ for (const auto& pt : quadRule) {
+
+ // get quadrature point
+ const Dune::FieldVector<ctype,dim>& quadPos = quad[pt].position();
+
+ // get integration factor
+ const ctype integrationElement = rGeom.integrationElement(quadPos);
+
+ // evaluate basis functions
+ dualCache->evaluateFunction(quadPos, values);
+
+ // compute values of function
+ T f_pos;
+ const GridFunction* gf = dynamic_cast<const GridFunction*>(&f_);
+ if (gf and gf->isDefinedOn(element))
+ gf->evaluateLocal(element, quadPos, f_pos);
+ else
+ f_.evaluate(geometry.global(quadPos), f_pos);
+
+ // and vector entries
+ for(size_t i=0; i<values.size(); ++i)
+ {
+ localVector[i].axpy(values[i]*quad[pt].weight()*integrationElement, f_pos);
+ }
+ }
+
+ for (size_t i=0; i<tFE.localBasis().size(); ++i) {
+ int idx = tBasis_.index(inside, i);
+ b[idx] += localB[i];
+ }
+}*/
int main(int argc, char *argv[]) { try {
Dune::MPIHelper::instance(argc, argv);
@@ -109,8 +230,8 @@ int main(int argc, char *argv[]) { try {
GlobalCoords upperRight0({2, 1});
buildGrid(lowerLeft0, upperRight0, n, grids[0]);
- GlobalCoords lowerLeft1({0.25, 1});
- GlobalCoords upperRight1({1.75, 2});
+ GlobalCoords lowerLeft1({0, 1});
+ GlobalCoords upperRight1({2, 2});
buildGrid(lowerLeft1, upperRight1, n, grids[1]);
// writing grids
@@ -139,27 +260,28 @@ int main(int argc, char *argv[]) { try {
CouplingPair coupling;
coupling.set(0, 1, upper, lower, 0.1, CouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
- double coveredArea_ = 0.8;
+ /* double coveredArea_ = 0.8;
CouplingType contactCoupling;
contactCoupling.setGrids(*grids[0], *grids[1]);
contactCoupling.setupContactPatch(*coupling.patch0(),*coupling.patch1());
contactCoupling.gridGlueBackend_ = coupling.backend();
contactCoupling.setCoveredArea(coveredArea_);
- contactCoupling.setup();
+ contactCoupling.setup();*/
// set nBodyAssembler
- /* std::vector<const GridType*> grids_ptr(grids.size());
+ using NBodyAssembler = Dune::Contact::NBodyAssembler<GridType, BlockVector>;
+ NBodyAssembler nBodyAssembler(grids.size(), 1);
+
+ std::vector<const GridType*> grids_ptr(grids.size());
for (size_t i=0; i<grids_ptr.size(); i++) {
grids_ptr[i] = grids[i].get();
}
- NBodyAssembler nBodyAssembler(grids.size(), 1);
nBodyAssembler.setGrids(grids_ptr);
-
nBodyAssembler.setCoupling(coupling, 0);
nBodyAssembler.assembleTransferOperator();
- nBodyAssembler.assembleObstacle(); */
+ nBodyAssembler.assembleObstacle();
// define basis
using Basis = P1NodalBasis<LevelView, double>;
@@ -177,142 +299,107 @@ int main(int argc, char *argv[]) { try {
printBasisDofLocation(basis1);
// set grid functions on coupling boundary
- BlockVector f0, f1;
+ std::vector<BlockVector> f(2);
BasisInterpolator<Basis, BlockVector, BoundaryFunction0<GlobalCoords>> interpolator0;
- interpolator0.interpolate(basis0, f0, BoundaryFunction0<GlobalCoords>());
+ interpolator0.interpolate(basis0, f[0], BoundaryFunction0<GlobalCoords>());
BasisInterpolator<Basis, BlockVector, BoundaryFunction1<GlobalCoords>> interpolator1;
- interpolator1.interpolate(basis1, f1, BoundaryFunction1<GlobalCoords>());
+ interpolator1.interpolate(basis1, f[1], BoundaryFunction1<GlobalCoords>());
- printBoundary(upper, basis0, f0, "f0: ");
- printBoundary(lower, basis1, f1, "f1: ");
- writeToVTK(basis0, f0, path, "body_0_level0");
- writeToVTK(basis1, f1, path, "body_1_level0");
+ BlockVector transformedF;
+ nBodyAssembler.nodalToTransformed(f, transformedF);
- // merged gridGlue coupling boundary
- auto& gridGlue = *contactCoupling.getGlue();
+ std::vector<Dune::BitSetVector<1>> boundaryVertices(2);
+ const auto& mortarCoupling = nBodyAssembler.getContactCouplings()[0];
+ const auto& nonmortarBoundary = mortarCoupling->nonmortarBoundary();
+ const auto& mortarBoundary = mortarCoupling->mortarBoundary();
- // make basis grid functions
- auto&& gridFunction0 = Functions::makeFunction(basis0, f0);
- auto&& gridFunction1 = Functions::makeFunction(basis1, f1);
+ nonmortarBoundary.getVertices(boundaryVertices[0]);
+ mortarBoundary.getVertices(boundaryVertices[1]);
- for (const auto& rIs : intersections(gridGlue)) {
- const auto& inside = rIs.inside();
- const auto& outside = rIs.outside();
-
- const auto& rGeometry = rIs.geometry();
- const auto& outGeometry = outside.geometry();
-
- for (size_t i=0; i<rGeometry.corners(); i++) {
- typename BasisGridFunction<Basis, BlockVector>::RangeType y;
- gridFunction1.evaluateLocal(outside, outGeometry.local(rGeometry.corner(i)), y);
- print(rGeometry.corner(i), "corner " + std::to_string(i));
- std::cout << "f1(corner) = " << y << std::endl;
-
- std::cout << std::endl;
- }
- std::cout << "---------"<< std::endl;
- // auto basis.
- // std::cout << containsInsideSubentity(rIs, 1, 1) << std::endl;
- }
-
- //Dune::GridGlue::GridGlueVtkWriter::writeIntersections<std::decltype(gridGlue), 0>(gridGlue, "");//, path + "gridglueintersections0.vtk"); //<std::decltype(gridGlue), 0>
-
- /*LevelBoundaryPatch boundaryWithMapping(grids[0]->levelGridView(0));
-
- const auto& indexSet0 = gridView0.indexSet();
- const auto& indexSet1 = gridView1.indexSet();
+ print(f[0], "f0: ");
+ print(boundaryVertices[0], "nonmortarBoundary: ");
+ print(f[1], "f1: ");
+ print(boundaryVertices[1], "mortarBoundary: ");
+ print(transformedF, "transformedF: ");
+ writeToVTK(basis0, f[0], path, "body_0_level0");
+ writeToVTK(basis1, f[1], path, "body_1_level0");
- // Get all fine grid boundary segments that are totally covered by the grid-glue segments
- typedef std::pair<int,int> Pair;
- std::map<Pair, c_type> coveredArea, fullArea;
+ print(mortarCoupling->mortarLagrangeMatrix(), "M: ");
+ print(mortarCoupling->nonmortarLagrangeMatrix(), "D: ");
- // initialize with area of boundary faces
- for (const auto& bIt : lower) {
- 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(gridGlue))
- coveredArea[Pair(indexSet0.index(rIs.inside()),rIs.indexInInside())] += rIs.geometry().volume();
+ std::vector<BlockVector> postprocessed(2);
+ nBodyAssembler.postprocess(transformedF, postprocessed);
- // add all fine grid faces that are totally covered by the contact mapping
- for (const auto& bIt : lower) {
- 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());
- }*/
+ print(postprocessed, "postprocessed: ");
- //writeBoundary(boundaryWithMapping, path + "relevantNonmortar");
+ /* using DualBasis = ;
+ DualBasis dualBasis;
+ BoundaryFunctionalAssembler<DualBasis> bAssembler(dualBasis, nonmortarBoundary);
+ ConstantFunction<GlobalCoords, Vector> oneFunction(1);
- /** \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());
+ BlockVector b;
+ L2FunctionalAssembler localAssembler(oneFunction);
+ bAssembler.assemble(localAssembler, b);
- /* printf("contact mapping could be built for %d of %d boundary segments.\n",
- boundaryWithMapping.numFaces(), lower.numFaces());
+ print(b, "b: ");*/
- LevelBoundaryPatch nonmortarBoundary = boundaryWithMapping;
- LevelBoundaryPatch mortarBoundary;
- mortarBoundary.setup(gridView1);
- for (const auto& rIs : intersections(gridGlue))
- if (nonmortarBoundary.contains(rIs.inside(),rIs.indexInInside()))
- mortarBoundary.addFace(rIs.outside(),rIs.indexInOutside());*/
+ /*
+ std::vector<BlockVector> g(2);
+ g[0].resize(f[0].size());
+ g[1].resize(f[1].size());
+ g[1][6] = 2;
+ BlockVector transformedG;
+ nBodyAssembler.nodalToTransformed(g, transformedG);
+ print(g[1], "g1: ");
+ print(transformedG, "transformedG: ");
+ // merged gridGlue coupling boundary
+ auto& gridGlue = *contactCoupling.getGlue();
+ // make basis grid functions
+ auto&& gridFunction0 = Functions::makeFunction(basis0, f[0]); */
+ /* for (const auto& bIs : intersections(upper)) {
+ const auto& inside = bIs.inside();
+ const auto& bGeometry = bIs.geometry();
- /** \todo Also restrict mortar indices and don't use the whole grid level. */
- /* Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));
+ for (size_t i=0; i<bGeometry.corners(); i++) {
- // 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)) {
+ typename BasisGridFunction<Basis, BlockVector>::RangeType y;
+ gridFunction1.evaluateLocal(outside, outGeometry.local(rGeometry.corner(i)), y);
+ print(rGeometry.corner(i), "corner " + std::to_string(i));
+ std::cout << "f1(corner) = " << y << std::endl;
+ }
+ }*/
- if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
- continue;
+ /*
+ auto&& gridFunction1 = Functions::makeFunction(basis1, f[1]);
+ for (const auto& rIs : intersections(gridGlue)) {
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;
+ const auto& rGeometry = rIs.geometry();
+ const auto& outGeometry = outside.geometry();
- for (int k=0; k<nTargetVertices; k++) {
- int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
- if (!mortarBoundary_.containsVertex(globalTargetIdx))
- continue;
+ for (size_t i=0; i<rGeometry.corners(); i++) {
+ typename BasisGridFunction<Basis, BlockVector>::RangeType y;
+ gridFunction1.evaluateLocal(outside, outGeometry.local(rGeometry.corner(i)), y);
+ print(rGeometry.corner(i), "corner " + std::to_string(i));
+ std::cout << "f1(corner) = " << y << std::endl;
- mortarIndices.add(localDomainIdx, globalTargetIdx);
- }
+ std::cout << std::endl;
}
- }
+ std::cout << "---------"<< std::endl;
+ }*/
- mortarIndices.exportIdx(mortarLagrangeMatrix_);
- // Clear it
- mortarLagrangeMatrix_ = 0; */
bool passed = true;
diff --git a/src/tests/nonoverlappingcouplingtest.cc b/src/tests/nonoverlappingcouplingtest.cc
new file mode 100644
index 00000000..83e5cee7
--- /dev/null
+++ b/src/tests/nonoverlappingcouplingtest.cc
@@ -0,0 +1,375 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+
+#include <array>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/geometrygrid.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+
+#include <dune/grid-glue/merging/contactmerge.hh>
+#include <dune/grid-glue/gridglue.hh>
+
+#include <dune/grid-glue/test/couplingtest.hh>
+#include <dune/grid-glue/test/communicationtest.hh>
+
+using namespace Dune;
+using namespace Dune::GridGlue;
+
+template<typename GridView>
+typename Dune::GridGlue::Codim1Extractor<GridView>::Predicate
+makeVerticalFacePredicate(double sliceCoord)
+{
+ using Element = typename GridView::Traits::template Codim<0>::Entity;
+ auto predicate = [sliceCoord](const Element& element, unsigned int face) -> bool {
+ const int dim = GridView::dimension;
+ const auto& refElement = Dune::ReferenceElements<double, dim>::general(element.type());
+
+ int numVertices = refElement.size(face, 1, dim);
+
+ for (int i=0; i<numVertices; i++)
+ if ( std::abs(element.geometry().corner(refElement.subEntity(face,1,i,dim))[0] - sliceCoord) > 1e-6 )
+ return false;
+
+ return true;
+ };
+ return predicate;
+}
+
+/** \brief trafo used for yaspgrids */
+template<int dim, typename ctype>
+class ShiftTrafo
+ : public AnalyticalCoordFunction< ctype, dim, dim, ShiftTrafo<dim,ctype> >
+{
+ double shift;
+public:
+ ShiftTrafo(double x) : shift(x) {};
+
+ //! evaluate method for global mapping
+ void evaluate ( const Dune::FieldVector<ctype, dim> &x, Dune::FieldVector<ctype, dim> &y ) const
+ {
+ y = x;
+ y[0] += shift;
+ }
+};
+
+template <int dim>
+void testMatchingCubeGrids()
+{
+
+ // ///////////////////////////////////////
+ // Make two cube grids
+ // ///////////////////////////////////////
+
+ using GridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+
+ std::array<int, dim> elements;
+ elements.fill(8);
+ FieldVector<double,dim> lower(0);
+ FieldVector<double,dim> upper(1);
+
+ GridType cubeGrid0(lower, upper, elements);
+
+ lower[0] += 1;
+ upper[0] += 1;
+
+ GridType cubeGrid1(lower, upper, elements);
+
+
+ // ////////////////////////////////////////
+ // Set up coupling at their interface
+ // ////////////////////////////////////////
+
+ typedef typename GridType::LevelGridView DomGridView;
+ typedef typename GridType::LevelGridView TarGridView;
+
+ typedef Codim1Extractor<DomGridView> DomExtractor;
+ typedef Codim1Extractor<TarGridView> TarExtractor;
+
+ const typename DomExtractor::Predicate domdesc = makeVerticalFacePredicate<DomGridView>(1);
+ const typename TarExtractor::Predicate tardesc = makeVerticalFacePredicate<TarGridView>(1);
+
+ auto domEx = std::make_shared<DomExtractor>(cubeGrid0.levelGridView(0), domdesc);
+ auto tarEx = std::make_shared<TarExtractor>(cubeGrid1.levelGridView(0), tardesc);
+
+ typedef Dune::GridGlue::GridGlue<DomExtractor,TarExtractor> GlueType;
+
+ // Testing with ContactMerge
+ typedef ContactMerge<dim,double> ContactMergeImpl;
+
+ auto contactMerger = std::make_shared<ContactMergeImpl>(0.01);
+ GlueType contactGlue(domEx, tarEx, contactMerger);
+
+ contactGlue.build();
+
+ std::cout << "Gluing successful, " << contactGlue.size() << " remote intersections found!" << std::endl;
+ assert(contactGlue.size() > 0);
+
+ // ///////////////////////////////////////////
+ // Test the coupling
+ // ///////////////////////////////////////////
+
+ testCoupling(contactGlue);
+ testCommunication(contactGlue);
+}
+
+
+template <int dim>
+void testNonMatchingCubeGrids()
+{
+
+ // ///////////////////////////////////////
+ // Make two cube grids
+ // ///////////////////////////////////////
+
+ using GridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+
+ std::array<int, dim> elements;
+ elements.fill(8);
+ FieldVector<double,dim> lower(0);
+ FieldVector<double,dim> upper(1);
+
+ GridType cubeGrid0(lower, upper, elements);
+
+ elements.fill(10);
+ lower[0] += 1;
+ upper[0] += 1;
+
+ GridType cubeGrid1(lower, upper, elements);
+
+
+ // ////////////////////////////////////////
+ // Set up coupling at their interface
+ // ////////////////////////////////////////
+
+ typedef typename GridType::LevelGridView DomGridView;
+ typedef typename GridType::LevelGridView TarGridView;
+
+ typedef Codim1Extractor<DomGridView> DomExtractor;
+ typedef Codim1Extractor<TarGridView> TarExtractor;
+
+ const typename DomExtractor::Predicate domdesc = makeVerticalFacePredicate<DomGridView>(1);
+ const typename TarExtractor::Predicate tardesc = makeVerticalFacePredicate<TarGridView>(1);
+
+ auto domEx = std::make_shared<DomExtractor>(cubeGrid0.levelGridView(0), domdesc);
+ auto tarEx = std::make_shared<TarExtractor>(cubeGrid1.levelGridView(0), tardesc);
+
+ typedef Dune::GridGlue::GridGlue<DomExtractor,TarExtractor> GlueType;
+
+ // Testing with ContactMerge
+ typedef ContactMerge<dim,double> ContactMergeImpl;
+
+ auto contactMerger = std::make_shared<ContactMergeImpl>(0.01);
+ GlueType contactGlue(domEx, tarEx, contactMerger);
+
+ contactGlue.build();
+
+ std::cout << "Gluing successful, " << contactGlue.size() << " remote intersections found!" << std::endl;
+ assert(contactGlue.size() > 0);
+
+ // ///////////////////////////////////////////
+ // Test the coupling
+ // ///////////////////////////////////////////
+
+ testCoupling(contactGlue);
+ testCommunication(contactGlue);
+}
+
+
+template<int dim, bool par>
+class MeshGenerator
+{
+ bool tar;
+public:
+ MeshGenerator(bool b) : tar(b) {};
+
+ using GridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+
+ std::shared_ptr<GridType> generate()
+ {
+ std::array<int, dim> elements;
+ elements.fill(8);
+ FieldVector<double,dim> lower(0);
+ FieldVector<double,dim> upper(1);
+
+ if (tar)
+ {
+ elements.fill(10);
+ lower[0] += 1;
+ upper[0] += 1;
+ }
+
+ return std::make_shared<GridType>(lower, upper, elements);
+ }
+};
+
+
+template<int dim>
+class MeshGenerator<dim, true>
+{
+ bool tar;
+public:
+ MeshGenerator(bool b) : tar(b) {};
+
+ typedef YaspGrid<dim> HostGridType;
+ typedef GeometryGrid<HostGridType, ShiftTrafo<dim,double> > GridType;
+
+ std::shared_ptr<GridType> generate()
+ {
+ std::array<int,dim> elements;
+ std::fill(elements.begin(), elements.end(), 8);
+ std::bitset<dim> periodic(0);
+ FieldVector<double,dim> size(1);
+ int overlap = 1;
+ double shift = 0.0;
+
+ if (tar)
+ {
+ std::fill(elements.begin(), elements.end(), 10);
+ shift = 1.0;
+ }
+
+ HostGridType * hostgridp = new HostGridType(
+ size, elements, periodic, overlap
+#if HAVE_MPI
+ , MPI_COMM_WORLD
+#endif
+ );
+ ShiftTrafo<dim,double> * trafop = new ShiftTrafo<dim,double>(shift);
+ return std::make_shared<GridType>(*hostgridp, *trafop);
+ }
+};
+
+
+template <int dim, class DomGen, class TarGen>
+void testParallelCubeGrids()
+{
+ // ///////////////////////////////////////
+ // Make two cube grids
+ // ///////////////////////////////////////
+
+ typedef typename DomGen::GridType GridType0;
+ typedef typename TarGen::GridType GridType1;
+
+ DomGen domGen(0);
+ TarGen tarGen(1);
+
+ double slice = 1.0;
+
+ std::shared_ptr<GridType0> cubeGrid0 = domGen.generate();
+ std::shared_ptr<GridType1> cubeGrid1 = tarGen.generate();
+
+ // ////////////////////////////////////////
+ // Set up Traits
+ // ////////////////////////////////////////
+
+ typedef typename GridType0::LevelGridView DomGridView;
+ typedef typename GridType1::LevelGridView TarGridView;
+
+ typedef Codim1Extractor<DomGridView> DomExtractor;
+ typedef Codim1Extractor<TarGridView> TarExtractor;
+
+ const typename DomExtractor::Predicate domdesc = makeVerticalFacePredicate<DomGridView>(slice);
+ const typename TarExtractor::Predicate tardesc = makeVerticalFacePredicate<TarGridView>(slice);
+
+ auto domEx = std::make_shared<DomExtractor>(cubeGrid0->levelGridView(0), domdesc);
+ auto tarEx = std::make_shared<TarExtractor>(cubeGrid1->levelGridView(0), tardesc);
+
+ // ////////////////////////////////////////
+ // Set up coupling at their interface
+ // ////////////////////////////////////////
+
+ typedef Dune::GridGlue::GridGlue<DomExtractor,TarExtractor> GlueType;
+
+ // Testing with ContactMerge
+ typedef ContactMerge<dim,double> ContactMergeImpl;
+
+ auto contactMerger = std::make_shared<ContactMergeImpl>(0.01);
+ GlueType contactGlue(domEx, tarEx, contactMerger);
+
+ contactGlue.build();
+
+ std::cout << "Gluing successful, " << contactGlue.size() << " remote intersections found!" << std::endl;
+ assert(contactGlue.size() > 0);
+
+ // ///////////////////////////////////////////
+ // Test the coupling
+ // ///////////////////////////////////////////
+
+ testCoupling(contactGlue);
+ testCommunication(contactGlue);
+}
+
+#if HAVE_MPI
+void eh( MPI_Comm *comm, int *err, ... )
+{
+ int len = 1024;
+ char error_txt[len];
+
+ MPI_Error_string(*err, error_txt, &len);
+ assert(len <= 1024);
+ DUNE_THROW(Dune::Exception, "MPI ERROR -- " << error_txt);
+}
+#endif // HAVE_MPI
+
+int main(int argc, char *argv[])
+{
+ Dune::MPIHelper::instance(argc, argv);
+ Dune::dinfo.attach(std::cout);
+
+#if HAVE_MPI
+ MPI_Errhandler errhandler;
+ MPI_Comm_create_errhandler(eh, &errhandler);
+ MPI_Comm_set_errhandler(MPI_COMM_WORLD, errhandler);
+#endif
+
+ // 2d Tests
+ typedef MeshGenerator<2,false> Seq;
+ typedef MeshGenerator<2,true> Par;
+
+ // Test two unit squares
+ std::cout << "==== 2D hybrid =============================================\n";
+ testMatchingCubeGrids<2>();
+ std::cout << "============================================================\n";
+ testNonMatchingCubeGrids<2>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<2,Seq,Seq>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<2,Par,Seq>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<2,Seq,Par>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<2,Par,Par>();
+ std::cout << "============================================================\n";
+
+ // 3d Tests
+#if ! HAVE_MPI
+ typedef MeshGenerator<3,false> Seq3d;
+ typedef MeshGenerator<3,true> Par3d;
+
+ // Test two unit cubes
+ std::cout << "==== 3D hybrid =============================================\n";
+ testMatchingCubeGrids<3>();
+ std::cout << "============================================================\n";
+ testNonMatchingCubeGrids<3>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<3,Seq3d,Seq3d>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<3,Par3d,Seq3d>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<3,Seq3d,Par3d>();
+ std::cout << "============================================================\n";
+ testParallelCubeGrids<3,Par3d,Par3d>();
+ std::cout << "============================================================\n";
+#endif // HAVE_MPI
+
+ return 0;
+}
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index d9618460..bfde1aeb 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -19,19 +19,23 @@ void IterationRegister::reset() {
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeStepper(
- Dune::ParameterTree const &parset,
- NBodyAssembler nBodyAssembler,
- GlobalFrictionContainer& globalFriction, Updaters ¤t,
- double relativeTime, double relativeTau,
- ExternalForces& externalForces,
- const std::vector<const ErrorNorm* >& errorNorms,
- std::function<bool(Updaters &, Updaters &)> mustRefine)
+ Dune::ParameterTree const &parset,
+ NBodyAssembler& nBodyAssembler,
+ GlobalFriction& globalFriction,
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+ Updaters ¤t,
+ double relativeTime,
+ double relativeTau,
+ ExternalForces& externalForces,
+ const std::vector<const ErrorNorm* >& errorNorms,
+ std::function<bool(Updaters &, Updaters &)> mustRefine)
: relativeTime_(relativeTime),
relativeTau_(relativeTau),
finalTime_(parset.get<double>("problem.finalTime")),
parset_(parset),
nBodyAssembler_(nBodyAssembler),
globalFriction_(globalFriction),
+ bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
current_(current),
R1_(),
externalForces_(externalForces),
@@ -100,7 +104,7 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(
Updaters const &oldUpdaters, double rTime, double rTau) {
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
- MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, globalFriction_,
+ MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, globalFriction_, bodywiseNonmortarBoundaries_,
newUpdatersAndCount.updaters, errorNorms_,
externalForces_)
.step(rTime, rTau);
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index 6e70c169..d0a014d0 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -28,15 +28,18 @@ class AdaptiveTimeStepper {
using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>;
- using GlobalFrictionContainer = typename MyCoupledTimeStepper::GlobalFrictionContainer;
+ using GlobalFriction = typename MyCoupledTimeStepper::GlobalFriction;
+ using BitVector = typename MyCoupledTimeStepper::BitVector;
using ExternalForces = typename MyCoupledTimeStepper::ExternalForces;
public:
AdaptiveTimeStepper(
Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
- GlobalFrictionContainer& globalFriction,
- Updaters ¤t, double relativeTime,
+ GlobalFriction& globalFriction,
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+ Updaters ¤t,
+ double relativeTime,
double relativeTau,
ExternalForces& externalForces,
const std::vector<const ErrorNorm* >& errorNorms,
@@ -55,7 +58,9 @@ class AdaptiveTimeStepper {
double finalTime_;
Dune::ParameterTree const &parset_;
NBodyAssembler& nBodyAssembler_;
- GlobalFrictionContainer& globalFriction_;
+ GlobalFriction& globalFriction_;
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
Updaters ¤t_;
UpdatersWithCount R1_;
ExternalForces& externalForces_;
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index a390e824..7908b085 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -5,29 +5,25 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
-#include <dune/common/function.hh>
-
#include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-#include <dune/tectonic/globalfriction.hh>
+#include "../spatial-solving/solverfactory_tmpl.cc"
+#include "../data-structures/levelcontactnetwork_tmpl.cc"
-#include "../data-structures/levelcontactnetwork.hh"
-#include "../spatial-solving/solverfactory.hh"
#include "rate/rateupdater.hh"
#include "state/stateupdater.hh"
#include "updaters.hh"
-/*
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
+using NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
+using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
-using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-template class AdaptiveTimeStepper<Factory, MyUpdaters, ErrorNorm>;
-*/
+using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
+
+template class AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index b9489572..0ee30c60 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -8,13 +8,16 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::CoupledTimeStepper(
double finalTime, Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
- GlobalFrictionContainer& globalFriction, Updaters updaters,
+ GlobalFriction& globalFriction,
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+ Updaters updaters,
const std::vector<const ErrorNorm* >& errorNorms,
ExternalForces& externalForces)
: finalTime_(finalTime),
parset_(parset),
nBodyAssembler_(nBodyAssembler),
globalFriction_(globalFriction),
+ bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
updaters_(updaters),
externalForces_(externalForces),
errorNorms_(errorNorms) {}
@@ -41,7 +44,7 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(double re
updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm> fixedPointIterator(
- parset_, nBodyAssembler_, globalFriction_, errorNorms_);
+ parset_, nBodyAssembler_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
velocityRHS, velocityIterate);
return iterations;
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 9d129d2d..bc2b1f67 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -14,14 +14,16 @@ class CoupledTimeStepper {
using Matrix = typename Factory::Matrix;
//using Nonlinearity = typename Factory::Nonlinearity;
public:
- using GlobalFrictionContainer = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::GlobalFrictionContainer;
+ using GlobalFriction = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::GlobalFriction;
+ using BitVector = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::BitVector;
using ExternalForces = std::vector<const std::function<void(double, Vector &)>*>;
public:
CoupledTimeStepper(double finalTime,
Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
- GlobalFrictionContainer& globalFriction,
+ GlobalFriction& globalFriction,
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
ExternalForces& externalForces);
@@ -31,7 +33,9 @@ class CoupledTimeStepper {
double finalTime_;
Dune::ParameterTree const &parset_;
NBodyAssembler& nBodyAssembler_;
- GlobalFrictionContainer& globalFriction_;
+ GlobalFriction& globalFriction_;
+ const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
Updaters updaters_;
ExternalForces& externalForces_;
const std::vector<const ErrorNorm* >& errorNorms_;
diff --git a/src/time-stepping/coupledtimestepper_tmpl.cc b/src/time-stepping/coupledtimestepper_tmpl.cc
index d5f6d1d1..c70a229f 100644
--- a/src/time-stepping/coupledtimestepper_tmpl.cc
+++ b/src/time-stepping/coupledtimestepper_tmpl.cc
@@ -5,24 +5,18 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
-#include <dune/common/function.hh>
-
#include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-#include <dune/tectonic/globalfriction.hh>
+#include "../spatial-solving/solverfactory_tmpl.cc"
+#include "../data-structures/levelcontactnetwork_tmpl.cc"
-#include "../data-structures/levelcontactnetwork.hh"
-#include "../spatial-solving/solverfactory.hh"
#include "rate/rateupdater.hh"
#include "state/stateupdater.hh"
#include "updaters.hh"
-/*
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
-
-using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+using NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
+using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
@@ -30,5 +24,5 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-template class CoupledTimeStepper<Factory, MyUpdaters, ErrorNorm>;
-*/
+
+template class CoupledTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
diff --git a/src/time-stepping/state.cc b/src/time-stepping/state.cc
index 02233a44..1e01e8a3 100644
--- a/src/time-stepping/state.cc
+++ b/src/time-stepping/state.cc
@@ -6,24 +6,43 @@
#include "state/ageinglawstateupdater.cc"
#include "state/sliplawstateupdater.cc"
-template <class ScalarVector, class Vector, class BoundaryPatchNodes>
+template <class ScalarVector, class Vector, class ContactCoupling, class FrictionCouplingPair>
auto initStateUpdater(
Config::stateModel model,
const std::vector<ScalarVector>& alpha_initial,
- const BoundaryPatchNodes& frictionNodes,
- const std::vector<double>& L,
- const std::vector<double>& V0)
+ const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+ const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings) // contains frictionInfo
-> std::shared_ptr<StateUpdater<ScalarVector, Vector>> {
+ assert(contactCouplings.size() == couplings.size());
+
+ auto stateUpdater = std::make_shared<StateUpdater<ScalarVector, Vector>>();
+
switch (model) {
case Config::AgeingLaw:
- return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionNodes, L, V0);
+ for (size_t i=0; i<couplings.size(); i++) {
+ const auto& coupling = couplings[i];
+ const auto nonmortarIdx = coupling->gridIdx_[0];
+
+ auto localUpdater = std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
+ alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+
+ stateUpdater->addLocalUpdater(localUpdater);
+ }
case Config::SlipLaw:
- return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionNodes, L, V0);
+ for (size_t i=0; i<couplings.size(); i++) {
+ const auto& coupling = couplings[i];
+ const auto nonmortarIdx = coupling->gridIdx_[0];
+
+ auto localUpdater = std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
+ alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+
+ stateUpdater->addLocalUpdater(localUpdater);
+ }
default:
assert(false);
+
+ return stateUpdater;
}
}
diff --git a/src/time-stepping/state.hh b/src/time-stepping/state.hh
index ca8352ad..81081ab8 100644
--- a/src/time-stepping/state.hh
+++ b/src/time-stepping/state.hh
@@ -3,19 +3,16 @@
#include <memory>
-#include <dune/common/bitsetvector.hh>
-
#include "../data-structures/enums.hh"
#include "state/ageinglawstateupdater.hh"
#include "state/sliplawstateupdater.hh"
#include "state/stateupdater.hh"
-template <class ScalarVector, class Vector, class BoundaryPatchNodes>
+template <class ScalarVector, class Vector, class ContactCoupling, class FrictionCouplingPair>
auto initStateUpdater(
Config::stateModel model,
const std::vector<ScalarVector>& alpha_initial,
- const BoundaryPatchNodes& frictionNodes,
- const std::vector<double>& L,
- const std::vector<double>& V0)
+ const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+ const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings) // contains frictionInfo
-> std::shared_ptr<StateUpdater<ScalarVector, Vector>>;
#endif
diff --git a/src/time-stepping/state/ageinglawstateupdater.cc b/src/time-stepping/state/ageinglawstateupdater.cc
index 94a7075c..6862029e 100644
--- a/src/time-stepping/state/ageinglawstateupdater.cc
+++ b/src/time-stepping/state/ageinglawstateupdater.cc
@@ -5,45 +5,37 @@
#include "../../utils/debugutils.hh"
template <class ScalarVector, class Vector>
-template <class BoundaryPatchNodes>
AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
- const std::vector<ScalarVector>& _alpha_initial,
- const BoundaryPatchNodes& _nodes,
- const std::vector<double>& _L,
- const std::vector<double>& _V0) :
- alpha(_alpha_initial),
- L(_L),
- V0(_V0) {
-
- size_t nBodies = alpha.size();
- nodes.resize(nBodies);
-
- for (size_t i=0; i<nBodies; i++) {
- auto& bodyNodes = nodes[i];
- const auto& _bcNodes = _nodes[i];
- if (_bcNodes.size() > 0) {
- bodyNodes.resize(_bcNodes[0]->size(), false);
- }
-
- for (size_t n=0; n<bodyNodes.size(); n++) {
- for (size_t bc=0; bc<_bcNodes.size(); bc++) {
- if (toBool((*_bcNodes[bc])[n])) {
- bodyNodes[n] = true;
- break;
- }
- }
- }
+ const ScalarVector& alpha_initial,
+ const BitVector& nodes,
+ const double L,
+ const double V0) :
+ nodes_(nodes),
+ L_(L),
+ V0_(V0) {
+
+ localToGlobal_.resize(nodes_.count());
+ alpha_.resize(localToGlobal_.size());
+
+ size_t localIdx = 0;
+ for (size_t i=0; i<nodes_.size(); i++) {
+ if (not toBool(nodes_[i]))
+ continue;
+
+ localToGlobal_[localIdx] = i;
+ alpha_[localIdx] = alpha_initial[i];
+ localIdx++;
}
}
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
- alpha_o = alpha;
+ alpha_o_ = alpha_;
}
template <class ScalarVector, class Vector>
-void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
- tau = _tau;
+void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double tau) {
+ tau_ = tau;
}
/*
@@ -61,36 +53,29 @@ auto liftSingularity(
}
template <class ScalarVector, class Vector>
-void AgeingLawStateUpdater<ScalarVector, Vector>::solve(const std::vector<Vector>& velocity_field) {
-
- for (size_t i = 0; i < alpha.size(); ++i) {
- const auto& velocity_fieldi = velocity_field[i];
- const auto& nodesi = nodes[i];
- auto& alphai = alpha[i];
- auto& alpha_oi = alpha_o[i];
-
- for (size_t j = 0; j < nodesi.size(); ++j) {
- if (not toBool(nodesi[j]))
- continue;
-
- double const V = velocity_fieldi[j].two_norm();
- double const mtoL = -tau / L[i];
- alphai[j] = std::log(std::exp(alpha_oi[j] + V * mtoL) +
- V0[i] * liftSingularity(mtoL, V));
- }
+void AgeingLawStateUpdater<ScalarVector, Vector>::solve(const Vector& velocity_field) {
+ for (size_t i=0; i<alpha_.size(); i++) {
+ double const V = velocity_field[i].two_norm();
+ double const mtoL = -tau_ / L_;
+ alpha_[i] = std::log(std::exp(alpha_o_[i] + V * mtoL) +
+ V0_ * liftSingularity(mtoL, V));
}
}
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
- std::vector<ScalarVector>& _alpha) {
+ ScalarVector& alpha) {
- _alpha = alpha;
+ assert(alpha.size() == nodes_.size());
+
+ for (size_t i=0; i<localToGlobal_.size(); i++) {
+ alpha[localToGlobal_[i]] = alpha_[i];
+ }
}
template <class ScalarVector, class Vector>
auto AgeingLawStateUpdater<ScalarVector, Vector>::clone() const
- -> std::shared_ptr<StateUpdater<ScalarVector, Vector>> {
+ -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> {
return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(*this);
}
diff --git a/src/time-stepping/state/ageinglawstateupdater.hh b/src/time-stepping/state/ageinglawstateupdater.hh
index 0c9ab450..58080a90 100644
--- a/src/time-stepping/state/ageinglawstateupdater.hh
+++ b/src/time-stepping/state/ageinglawstateupdater.hh
@@ -1,31 +1,37 @@
#ifndef SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
#define SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
+#include <dune/common/bitsetvector.hh>
+
#include "stateupdater.hh"
template <class ScalarVector, class Vector>
-class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
+class AgeingLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+ using BitVector = Dune::BitSetVector<1>;
+
public:
- template <class BoundaryPatchNodes>
AgeingLawStateUpdater(
- const std::vector<ScalarVector>& _alpha_initial,
- const BoundaryPatchNodes& _nodes,
- const std::vector<double>& _L,
- const std::vector<double>& _V0);
+ const ScalarVector& alpha_initial,
+ const BitVector& nodes,
+ const double L,
+ const double V0);
void nextTimeStep() override;
- void setup(double _tau) override;
- void solve(const std::vector<Vector>& velocity_field) override;
- void extractAlpha(std::vector<ScalarVector> &) override;
+ void setup(double tau) override;
+ void solve(const Vector& velocity_field) override;
+ void extractAlpha(ScalarVector&) override;
- auto clone() const -> std::shared_ptr<StateUpdater<ScalarVector, Vector>> override;
+ auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
private:
- std::vector<ScalarVector> alpha_o;
- std::vector<ScalarVector> alpha;
- std::vector<Dune::BitSetVector<1>> nodes;
- const std::vector<double>& L;
- const std::vector<double>& V0;
- double tau;
+ std::vector<int> localToGlobal_;
+
+ ScalarVector alpha_o_;
+ ScalarVector alpha_;
+ const BitVector& nodes_;
+ const double L_;
+ const double V0_;
+ double tau_;
};
#endif
diff --git a/src/time-stepping/state/sliplawstateupdater.cc b/src/time-stepping/state/sliplawstateupdater.cc
index 51ae920b..3b06f055 100644
--- a/src/time-stepping/state/sliplawstateupdater.cc
+++ b/src/time-stepping/state/sliplawstateupdater.cc
@@ -4,78 +4,64 @@
#include "../../utils/tobool.hh"
template <class ScalarVector, class Vector>
-template <class BoundaryPatchNodes>
SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
- const std::vector<ScalarVector>& _alpha_initial,
- const BoundaryPatchNodes& _nodes,
- const std::vector<double>& _L,
- const std::vector<double>& _V0) :
- alpha(_alpha_initial),
- L(_L),
- V0(_V0) {
+ const ScalarVector& alpha_initial,
+ const BitVector& nodes,
+ const double L,
+ const double V0) :
+ nodes_(nodes),
+ L_(L),
+ V0_(V0) {
- size_t nBodies = alpha.size();
- nodes.resize(nBodies);
+ localToGlobal_.resize(nodes_.count());
+ alpha_.resize(localToGlobal_.size());
- for (size_t i=0; i<nBodies; i++) {
- auto& bodyNodes = nodes[i];
- auto& _bodyNodes = _nodes[i];
- if (_bodyNodes.size() > 0) {
- bodyNodes.resize(_bodyNodes[0]->size(), false);
- }
+ size_t localIdx = 0;
+ for (size_t i=0; i<nodes_.size(); i++) {
+ if (not toBool(nodes_[i]))
+ continue;
- for (size_t n=0; n<bodyNodes.size(); n++) {
- for (size_t bc=0; bc<_bodyNodes.size(); bc++) {
- if (toBool((*_bodyNodes[bc])[n])) {
- bodyNodes[n] = true;
- break;
- }
- }
- }
+ localToGlobal_[localIdx] = i;
+ alpha_[localIdx] = alpha_initial[i];
+ localIdx++;
}
}
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
- alpha_o = alpha;
+ alpha_o_ = alpha_;
}
template <class ScalarVector, class Vector>
-void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
- tau = _tau;
+void SlipLawStateUpdater<ScalarVector, Vector>::setup(double tau) {
+ tau_ = tau;
}
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::solve(
- const std::vector<Vector>& velocity_field) {
+ const Vector& velocity_field) {
- for (size_t i = 0; i < alpha.size(); ++i) {
- const auto& velocity_fieldi = velocity_field[i];
- const auto& nodesi = nodes[i];
- auto& alphai = alpha[i];
- auto& alpha_oi = alpha_o[i];
-
- for (size_t j = 0; j < nodesi.size(); ++j) {
- if (not toBool(nodesi[j]))
- continue;
-
- double const V = velocity_fieldi[j].two_norm();
- double const mtoL = -tau * V / L[i];
- alphai[j] = (V <= 0) ? alpha_oi[j] : std::expm1(mtoL) * std::log(V / V0[i]) +
- alpha_oi[j] * std::exp(mtoL);
- }
+ for (size_t i=0; i<alpha_.size(); i++) {
+ double const V = velocity_field[i].two_norm();
+ double const mtVoL = -tau_ * V / L_;
+ alpha_[i] = (V <= 0) ? alpha_o_[i] : std::expm1(mtVoL) * std::log(V / V0_) +
+ alpha_o_[i] * std::exp(mtVoL);
}
}
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
- std::vector<ScalarVector>& _alpha) {
+ ScalarVector& alpha) {
+
+ assert(alpha.size() == nodes_.size());
- _alpha = alpha;
+ for (size_t i=0; i<localToGlobal_.size(); i++) {
+ alpha[localToGlobal_[i]] = alpha_[i];
+ }
}
template <class ScalarVector, class Vector>
auto SlipLawStateUpdater<ScalarVector, Vector>::clone() const
--> std::shared_ptr<StateUpdater<ScalarVector, Vector>> {
+-> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> {
return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(*this);
-}
+}
diff --git a/src/time-stepping/state/sliplawstateupdater.hh b/src/time-stepping/state/sliplawstateupdater.hh
index 0b3f7ef3..478a8b62 100644
--- a/src/time-stepping/state/sliplawstateupdater.hh
+++ b/src/time-stepping/state/sliplawstateupdater.hh
@@ -1,32 +1,38 @@
#ifndef SRC_TIME_STEPPING_STATE_SLIPLAWSTATEUPDATER_HH
#define SRC_TIME_STEPPING_STATE_SLIPLAWSTATEUPDATER_HH
+#include <dune/common/bitsetvector.hh>
+
#include "stateupdater.hh"
template <class ScalarVector, class Vector>
-class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
+class SlipLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+ using BitVector = Dune::BitSetVector<1>;
+
public:
- template <class BoundaryPatchNodes>
SlipLawStateUpdater(
- const std::vector<ScalarVector>& _alpha_initial,
- const BoundaryPatchNodes& _nodes,
- const std::vector<double>& _L,
- const std::vector<double>& _V0);
+ const ScalarVector& alpha_initial,
+ const BitVector& nodes,
+ const double L,
+ const double V0);
- void nextTimeStep() override;
- void setup(double _tau) override;
- void solve(const std::vector<Vector>& velocity_field) override;
- void extractAlpha(std::vector<ScalarVector> &) override;
+ void nextTimeStep() override;
+ void setup(double tau) override;
+ void solve(const Vector& velocity_field) override;
+ void extractAlpha(ScalarVector&) override;
- auto clone() const -> std::shared_ptr<StateUpdater<ScalarVector, Vector>> override;
+ auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
private:
- std::vector<ScalarVector> alpha_o;
- std::vector<ScalarVector> alpha;
- std::vector<Dune::BitSetVector<1>> nodes;
- const std::vector<double>& L;
- const std::vector<double>& V0;
- double tau;
+ std::vector<int> localToGlobal_;
+
+ ScalarVector alpha_o_;
+ ScalarVector alpha_;
+ const BitVector& nodes_;
+ const double L_;
+ const double V0_;
+ double tau_;
};
#endif
diff --git a/src/time-stepping/state/stateupdater.hh b/src/time-stepping/state/stateupdater.hh
index f7dd1b1b..ed7a5241 100644
--- a/src/time-stepping/state/stateupdater.hh
+++ b/src/time-stepping/state/stateupdater.hh
@@ -1,16 +1,74 @@
#ifndef SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
#define SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
-template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
+#include <memory>
+#include <vector>
+
+// state updater for each coupling
+template <class ScalarVectorTEMPLATE, class Vector> class LocalStateUpdater {
public:
using ScalarVector = ScalarVectorTEMPLATE;
+ void setBodyIndex(const size_t bodyIdx) {
+ bodyIdx_ = bodyIdx;
+ }
+
+ auto bodyIndex() const {
+ return bodyIdx_;
+ }
+
void virtual nextTimeStep() = 0;
void virtual setup(double _tau) = 0;
- void virtual solve(const std::vector<Vector>& velocity_field) = 0;
- void virtual extractAlpha(std::vector<ScalarVector>& alpha) = 0;
+ void virtual solve(const Vector& velocity_field) = 0;
+ void virtual extractAlpha(ScalarVector& alpha) = 0;
+
+ std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
- std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
+private:
+ size_t bodyIdx_;
};
+
+template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
+public:
+ using ScalarVector = ScalarVectorTEMPLATE;
+ using LocalStateUpdater = LocalStateUpdater<ScalarVector, Vector>;
+
+ void addLocalUpdater(std::shared_ptr<LocalStateUpdater> localStateUpdater) {
+ localStateUpdaters_.emplace_back(localStateUpdater);
+ }
+
+ void nextTimeStep() {
+ for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+ localStateUpdaters_[i]->nextTimeStep();
+ }
+ }
+
+ void setup(double tau) {
+ for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+ localStateUpdaters_[i]->setup(tau);
+ }
+ }
+
+ void solve(const std::vector<Vector>& velocity_field) {
+ for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+ auto& localStateUpdater = localStateUpdaters_[i];
+ localStateUpdater->solve(velocity_field[localStateUpdater->bodyIndex()]);
+ }
+ }
+
+ void extractAlpha(std::vector<ScalarVector>& alpha) {
+ for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+ auto& localStateUpdater = localStateUpdaters_[i];
+ localStateUpdater->extractAlpha(alpha[localStateUpdater->bodyIndex()]);
+ }
+ }
+
+ std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const {
+ return std::make_shared<StateUpdater<ScalarVector, Vector>>(*this);
+ }
+
+private:
+ std::vector<std::shared_ptr<LocalStateUpdater>> localStateUpdaters_;
+};
#endif
diff --git a/src/time-stepping/state_tmpl.cc b/src/time-stepping/state_tmpl.cc
index 42fa0b33..39fd6bd4 100644
--- a/src/time-stepping/state_tmpl.cc
+++ b/src/time-stepping/state_tmpl.cc
@@ -1,19 +1,20 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
#include "../explicitvectors.hh"
#include "../explicitgrid.hh"
-#include "../data-structures/levelcontactnetwork.hh"
+#include <dune/common/promotiontraits.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
+#include "../frictioncouplingpair.hh"
+
+using field_type = typename Dune::PromotionTraits<typename Vector::field_type,
+ typename DeformedGrid::ctype>::PromotedType;
-using BoundaryPatchNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryPatchNodes;
+using MyContactCoupling = Dune::Contact::DualMortarCoupling<field_type, DeformedGrid>;
+using MyFrictionCouplingPair = FrictionCouplingPair<DeformedGrid, LocalVector, field_type>;
-template
-auto initStateUpdater<ScalarVector, Vector, BoundaryPatchNodes>(
- Config::stateModel model,
- const std::vector<ScalarVector>& alpha_initial,
- const BoundaryPatchNodes& frictionalNodes,
- const std::vector<double>& L,
- const std::vector<double>& V0)
--> std::shared_ptr<StateUpdater<ScalarVector, Vector>>;
+template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+initStateUpdater<ScalarVector, Vector, MyContactCoupling, MyFrictionCouplingPair>(
+ Config::stateModel model,
+ const std::vector<ScalarVector>& alpha_initial,
+ const std::vector<std::shared_ptr<MyContactCoupling>>& contactCouplings,
+ const std::vector<std::shared_ptr<MyFrictionCouplingPair>>& couplings);
diff --git a/src/utils/almostequal.hh b/src/utils/almostequal.hh
index 448c95fe..d75222d7 100644
--- a/src/utils/almostequal.hh
+++ b/src/utils/almostequal.hh
@@ -1,6 +1,8 @@
#ifndef SRC_ALMOSTEQUAL_HH
#define SRC_ALMOSTEQUAL_HH
+#include <type_traits>
+#include <limits>
#include <math.h>
template <typename ctype>
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index c81de410..cb289f16 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -23,6 +23,7 @@
for (size_t j=0; j<s; j++) {
std::cout << x[i][j] << " ";
}
+ std::cout << std::endl;
}
std::cout << std::endl << std::endl;
}
--
GitLab