diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fbd99d9fe11f61a3ea84653dbbf8cb71bf7ee551..d86dde2a5b076bb28e0578c31c69df606c69ca8f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory("tests")
+
set(SW_SOURCE_FILES
assemblers.cc
data-structures/body.cc
@@ -24,10 +26,10 @@ set(SW_SOURCE_FILES
set(MSW_SOURCE_FILES
assemblers.cc
data-structures/body.cc
- data-structures/contactnetwork.cc
+ #data-structures/contactnetwork.cc
data-structures/enumparser.cc
data-structures/levelcontactnetwork.cc
- factories/cantorfactory.cc
+ #factories/cantorfactory.cc
factories/stackedblocksfactory.cc
io/vtk.cc
io/hdf5/frictionalboundary-writer.cc
@@ -36,8 +38,11 @@ set(MSW_SOURCE_FILES
io/hdf5/restart-io.cc
io/hdf5/surface-writer.cc
io/hdf5/time-writer.cc
- multi-body-problem-data/cuboidgeometry.cc
- multi-body-problem-data/mygrids.cc
+ #multi-body-problem-data/grid/cube.cc
+ #multi-body-problem-data/grid/cubefaces.cc
+ multi-body-problem-data/grid/cuboidgeometry.cc
+ 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
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 1097625fbd3055c9b7d68720d93b5fbc3952c2c1..d19b8e32f1f77790ad491479fadb4f97dba28634 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -24,8 +24,8 @@
#include "assemblers.hh"
template <class GridView, int dimension>
-MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
- : cellBasis(_gridView),
+MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView) :
+ cellBasis(_gridView),
vertexBasis(_gridView),
gridView(_gridView),
cellAssembler(cellBasis, cellBasis),
@@ -33,8 +33,9 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
- BoundaryPatch<GridView> const &frictionalBoundary,
- ScalarMatrix &frictionalBoundaryMass) const {
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ ScalarMatrix &frictionalBoundaryMass) const {
+
BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
frictionalBoundaryMassAssembler(frictionalBoundary);
@@ -44,9 +45,9 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleMass(
- Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
- densityFunction,
- Matrix &M) const {
+ Dune::VirtualFunction<LocalVector, LocalScalarVector> const & densityFunction,
+ Matrix &M) const {
+
// NOTE: We treat the weight as a constant function
QuadratureRuleKey quadKey(dimension, 0);
@@ -58,8 +59,11 @@ void MyAssembler<GridView, dimension>::assembleMass(
}
template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
- Matrix &A) const {
+void MyAssembler<GridView, dimension>::assembleElasticity(
+ double E,
+ double nu,
+ Matrix &A) const {
+
StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
localStiffness(E, nu);
vertexAssembler.assembleOperator(localStiffness, A);
@@ -67,9 +71,10 @@ void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleViscosity(
- Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
- Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
- Matrix &C) const {
+ Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
+ Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
+ Matrix &C) const {
+
// NOTE: We treat the weights as constant functions
QuadratureRuleKey shearViscosityKey(dimension, 0);
QuadratureRuleKey bulkViscosityKey(dimension, 0);
@@ -83,8 +88,9 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleBodyForce(
- Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
- Vector &f) const {
+ Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+ Vector &f) const {
+
L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
gravityFunctionalAssembler(gravityField);
vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
@@ -92,9 +98,11 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNeumann(
- BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
- Dune::VirtualFunction<double, double> const &neumann,
- double relativeTime) const {
+ BoundaryPatch<GridView> const &neumannBoundary,
+ Vector &f,
+ Dune::VirtualFunction<double, double> const &neumann,
+ double relativeTime) const {
+
LocalVector localNeumann(0);
neumann.evaluate(relativeTime, localNeumann[0]);
NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
@@ -106,9 +114,11 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
- BoundaryPatch<GridView> const &frictionalBoundary,
- ScalarVector &weightedNormalStress, double youngModulus,
- double poissonRatio, Vector const &displacement) const {
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ ScalarVector &weightedNormalStress,
+ double youngModulus,
+ double poissonRatio,
+ Vector const &displacement) const {
BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
displacement);
@@ -138,11 +148,12 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
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>> {
+ 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;
{
@@ -169,8 +180,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleVonMisesStress(
- double youngModulus, double poissonRatio, Vector const &u,
- ScalarVector &stress) const {
+ double youngModulus,
+ double poissonRatio,
+ Vector const &u,
+ ScalarVector &stress) const {
+
auto const gridDisplacement =
std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
vertexBasis, u);
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 195370d454400ba28d9a20ca9dc5630db916c124..a32870a9d25159d8d00420323e002d4f75b0607c 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -18,73 +18,83 @@
#include "data-structures/enums.hh"
-template <class GridView, int dimension> class MyAssembler {
+template <class GridView, int dimension>
+class MyAssembler {
public:
- using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
- using Matrix =
+ using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+ using Matrix =
Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
- using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
- using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+ using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+ using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
- using CellBasis = P0Basis<GridView, double>;
- using VertexBasis = P1NodalBasis<GridView, double>;
+ using CellBasis = P0Basis<GridView, double>;
+ using VertexBasis = P1NodalBasis<GridView, double>;
- CellBasis const cellBasis;
- VertexBasis const vertexBasis;
- GridView const &gridView;
+ CellBasis const cellBasis;
+ VertexBasis const vertexBasis;
+ GridView const &gridView;
private:
- using Grid = typename GridView::Grid;
- using LocalVector = typename Vector::block_type;
- using LocalScalarVector = typename ScalarVector::block_type;
+ using Grid = typename GridView::Grid;
+ using LocalVector = typename Vector::block_type;
+ using LocalScalarVector = typename ScalarVector::block_type;
- using LocalCellBasis = typename CellBasis::LocalFiniteElement;
- using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
+ using LocalCellBasis = typename CellBasis::LocalFiniteElement;
+ using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
- Assembler<CellBasis, CellBasis> cellAssembler;
- Assembler<VertexBasis, VertexBasis> vertexAssembler;
+ Assembler<CellBasis, CellBasis> cellAssembler;
+ Assembler<VertexBasis, VertexBasis> vertexAssembler;
public:
- MyAssembler(GridView const &gridView);
-
- void assembleFrictionalBoundaryMass(
- BoundaryPatch<GridView> const &frictionalBoundary,
- ScalarMatrix &frictionalBoundaryMass) const;
-
- void assembleMass(Dune::VirtualFunction<
- LocalVector, LocalScalarVector> const &densityFunction,
- Matrix &M) const;
-
- void assembleElasticity(double E, double nu, Matrix &A) const;
-
- void assembleViscosity(
- Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
- shearViscosity,
- Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
- bulkViscosity,
- Matrix &C) const;
-
- void assembleBodyForce(
- Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
- Vector &f) const;
-
- void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
- Vector &f,
- Dune::VirtualFunction<double, double> const &neumann,
- double relativeTime) const;
-
- void assembleWeightedNormalStress(
- BoundaryPatch<GridView> const &frictionalBoundary,
- ScalarVector &weightedNormalStress, double youngModulus,
- double poissonRatio, Vector const &displacement) const;
-
- std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
- Config::FrictionModel frictionModel,
- BoundaryPatch<GridView> const &frictionalBoundary,
- GlobalFrictionData<dimension> const &frictionInfo,
- ScalarVector const &weightedNormalStress) const;
-
- void assembleVonMisesStress(double youngModulus, double poissonRatio,
- Vector const &u, ScalarVector &stress) const;
+ MyAssembler(GridView const &gridView);
+
+ void assembleFrictionalBoundaryMass(
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ ScalarMatrix &frictionalBoundaryMass) const;
+
+ void assembleMass(
+ Dune::VirtualFunction<
+ LocalVector, LocalScalarVector> const &densityFunction,
+ Matrix &M) const;
+
+ void assembleElasticity(
+ double E,
+ double nu,
+ Matrix &A) const;
+
+ void assembleViscosity(
+ Dune::VirtualFunction<LocalVector, LocalScalarVector> const & shearViscosity,
+ Dune::VirtualFunction<LocalVector, LocalScalarVector> const & bulkViscosity,
+ Matrix &C) const;
+
+ void assembleBodyForce(
+ Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+ Vector &f) const;
+
+ void assembleNeumann(
+ BoundaryPatch<GridView> const &neumannBoundary,
+ Vector &f,
+ Dune::VirtualFunction<double, double> const &neumann,
+ double relativeTime) const;
+
+ void assembleWeightedNormalStress(
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ ScalarVector &weightedNormalStress,
+ double youngModulus,
+ double poissonRatio,
+ Vector const &displacement) const;
+
+ auto assembleFrictionNonlinearity(
+ Config::FrictionModel frictionModel,
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ GlobalFrictionData<dimension> const &frictionInfo,
+ ScalarVector const &weightedNormalStress) const
+ -> std::shared_ptr<GlobalFriction<Matrix, Vector>>;
+
+ void assembleVonMisesStress(
+ double youngModulus,
+ double poissonRatio,
+ Vector const &u,
+ ScalarVector &stress) const;
};
#endif
diff --git a/src/data-structures/body.cc b/src/data-structures/body.cc
index f110eec13b143b1b619a35c5ac4f248cef133785..b72f396b28523faf0e48e46fe63510e2edb835f2 100644
--- a/src/data-structures/body.cc
+++ b/src/data-structures/body.cc
@@ -6,7 +6,9 @@
#include "body.hh"
template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-Body<GridTEMPLATE, VectorTEMPLATE, dims>::Body(const std::shared_ptr<BodyData<dims>>& bodyData, const std::shared_ptr<GridType>& grid) :
+Body<GridTEMPLATE, VectorTEMPLATE, dims>::Body(
+ const std::shared_ptr<BodyData<dims>>& bodyData,
+ const std::shared_ptr<GridType>& grid) :
bodyData_(bodyData),
grid_(grid),
deformation_(std::make_shared<DeformationFunction>(grid_->leafGridView())),
@@ -22,6 +24,23 @@ Body<GridTEMPLATE, VectorTEMPLATE, dims>::Body(const std::shared_ptr<BodyData<di
for (size_t i=0; i<levelBoundaryConditions_.size(); i++) {
levelBoundaryConditions_[i].resize(0);
}
+
+ externalForce_ = [&](double relativeTime, VectorType &ell) {
+ // Problem formulation: right-hand side
+ std::vector<std::shared_ptr<LeafBoundaryCondition>> leafNeumannConditions;
+ leafBoundaryConditions("neumann", leafNeumannConditions);
+ ell.resize(gravityFunctional_.size());
+ ell = gravityFunctional_;
+
+ for (size_t i=0; i<leafNeumannConditions.size(); i++) {
+ const auto& leafNeumannCondition = leafNeumannConditions[i];
+
+ VectorType b;
+ assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(), b, *leafNeumannCondition->boundaryFunction(),
+ relativeTime);
+ ell += b;
+ }
+ };
}
@@ -51,8 +70,11 @@ void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assemble() {
}
template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembleFriction(const Config::FrictionModel& frictionModel, const ScalarVector& weightedNormalStress) {
- globalFriction_.resize(frictionBoundaryConditions_.size());
+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];
@@ -60,21 +82,241 @@ void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembleFriction(const Config::Fr
}
}
+// setter and getter
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::data() const
+-> const std::shared_ptr<BodyData<dims>>& {
+
+ return bodyData_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::setDeformation(const VectorType& def) {
+ deformation_->setDeformation(def);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::deformation() const
+-> DeformationFunction& {
+
+ return *deformation_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::deformedGrid() const
+-> const DeformedGridType* {
+
+ return deformedGrid_.get();
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafView() const
+-> const DeformedLeafGridView& {
+
+ return leafView_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafVertexCount() const
+-> int {
+
+ return leafVertexCount_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembler() const
+-> const std::shared_ptr<Assembler>& {
+
+ return assembler_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::matrices() const
+-> const Matrices<Matrix,1>& {
+
+ return matrices_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::externalForce() const
+-> const ExternalForce& {
+
+ 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);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+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]);
+ }
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafBoundaryConditions() const
+-> const std::vector<std::shared_ptr<LeafBoundaryCondition>>& {
+
+ return leafBoundaryConditions_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(
+ std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition,
+ int level) {
+
+ levelBoundaryConditions_[level].push_back(levelBoundaryCondition);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions(
+ const std::string& tag,
+ std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const {
+
+ selectedConditions.resize(levelBoundaryConditions_.size());
+ for (size_t level=0; level<levelBoundaryConditions_.size(); level++) {
+ const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
+ for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
+ if (levelBoundaryConditions[i]->tag() == tag)
+ selectedConditions[level].push_back(levelBoundaryConditions[i]);
+ }
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions(
+ const std::string& tag,
+ int level,
+ std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const {
+
+ selectedConditions.resize(0);
+ const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
+ for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
+ if (levelBoundaryConditions[i]->tag() == tag)
+ selectedConditions.push_back(levelBoundaryConditions[i]);
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions(int level) const
+-> const std::vector<std::shared_ptr<LevelBoundaryCondition>>& {
+
+ return levelBoundaryConditions_[level];
+}
+
template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::externalForce(double relativeTime, VectorType& ell) const {
- // Problem formulation: right-hand side
- std::vector<std::shared_ptr<LeafBoundaryCondition>> leafNeumannConditions;
- leafBoundaryConditions("neumann", leafNeumannConditions);
- ell.resize(gravityFunctional_.size());
- ell = gravityFunctional_;
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions() const
+-> const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& {
+
+ return levelBoundaryConditions_;
+}
+
+
+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,
+ BoundaryPatchNodes& nodes) const {
+
+ std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+ this->leafBoundaryConditions(tag, boundaryConditions);
+
+ nodes.resize(boundaryConditions.size());
+
+ for (size_t i=0; i<nodes.size(); i++) {
+ nodes[i] = boundaryConditions[i]->boundaryPatch()->getVertices();
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryNodes(
+ const std::string& tag,
+ BoundaryNodes& nodes) const {
+
+ std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+ this->leafBoundaryConditions(tag, boundaryConditions);
+
+ nodes.resize(boundaryConditions.size());
+
+ for (size_t i=0; i<nodes.size(); i++) {
+ nodes[i] = boundaryConditions[i]->boundaryNodes().get();
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryFunctions(
+ const std::string& tag,
+ BoundaryFunctions& functions) const {
+
+ std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+ this->leafBoundaryConditions(tag, boundaryConditions);
+
+ functions.resize(boundaryConditions.size());
+
+ for (size_t i=0; i<functions.size(); i++) {
+ functions[i] = boundaryConditions[i]->boundaryFunction().get();
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryPatches(
+ const std::string& tag,
+ BoundaryPatches& patches) const {
+
+ std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+ this->leafBoundaryConditions(tag, boundaryConditions);
- for (size_t i=0; i<leafNeumannConditions.size(); i++) {
- const auto& leafNeumannCondition = leafNeumannConditions[i];
+ patches.resize(boundaryConditions.size());
- VectorType b;
- assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(), b, *leafNeumannCondition->boundaryFunction(),
- relativeTime);
- ell += b;
+ for (size_t i=0; i<patches.size(); i++) {
+ patches[i] = boundaryConditions[i]->boundaryPatch().get();
}
}
diff --git a/src/data-structures/body.hh b/src/data-structures/body.hh
index 5359fbd096bb701967f70b319a0de6bc395c051a..6a254e4ccb1e92e75216bf3db6b8166f4a3a2490 100644
--- a/src/data-structures/body.hh
+++ b/src/data-structures/body.hh
@@ -1,5 +1,5 @@
-#ifndef SRC_BODY_HH
-#define SRC_BODY_HH
+#ifndef SRC_DATA_STRUCTURES_BODY_HH
+#define SRC_DATA_STRUCTURES_BODY_HH
#include <functional>
@@ -29,152 +29,136 @@
#include "../assemblers.hh"
#include "../boundarycondition.hh"
+#include "globalfrictioncontainer.hh"
#include "enums.hh"
#include "matrices.hh"
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims> class Body {
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+class Body {
public:
- using GridType = GridTEMPLATE;
- using VectorType = VectorTEMPLATE;
+ using GridType = GridTEMPLATE;
+ using VectorType = VectorTEMPLATE;
- using DeformationFunction = DeformationFunction<typename GridType::LeafGridView, VectorType>;
- using DeformedGridType = Dune::GeometryGrid<GridType, DeformationFunction>;
- using DeformedLeafGridView = typename DeformedGridType::LeafGridView;
- using DeformedLevelGridView = typename DeformedGridType::LevelGridView;
+ using DeformationFunction = DeformationFunction<typename GridType::LeafGridView, VectorType>;
+ using DeformedGridType = Dune::GeometryGrid<GridType, DeformationFunction>;
+ using DeformedLeafGridView = typename DeformedGridType::LeafGridView;
+ using DeformedLevelGridView = typename DeformedGridType::LevelGridView;
- using Assembler = MyAssembler<DeformedLeafGridView, dims>;
- using Matrix = typename Assembler::Matrix;
- using LocalVector = typename VectorType::block_type;
+ using Assembler = MyAssembler<DeformedLeafGridView, dims>;
+ using Matrix = typename Assembler::Matrix;
+ using LocalVector = typename VectorType::block_type;
- using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
- using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+ using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+ using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
- using Function = Dune::VirtualFunction<double, double>;
+ 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 LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
+ using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
+ using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
+ using LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
+
+ using BoundaryFunctions = std::vector<const Function* >;
+ using BoundaryNodes = std::vector<const Dune::BitSetVector<dims>* >;
+ 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:
- std::shared_ptr<BodyData<dims>> bodyData_;
+ std::shared_ptr<Body<GridType, VectorType, dims>> parent_;
+ std::vector<std::shared_ptr<Body<GridType, VectorType, dims>>> children_;
+
+ std::shared_ptr<BodyData<dims>> bodyData_;
+
+ std::shared_ptr<GridType> grid_;
+ std::shared_ptr<DeformationFunction> deformation_;
+ std::shared_ptr<DeformedGridType> deformedGrid_;
- std::shared_ptr<GridType> grid_;
- std::shared_ptr<DeformationFunction> deformation_;
- std::shared_ptr<DeformedGridType> deformedGrid_;
+ DeformedLeafGridView leafView_;
+ int leafVertexCount_;
- DeformedLeafGridView leafView_;
- int leafVertexCount_;
+ std::shared_ptr<Assembler> assembler_;
+ Matrices<Matrix,1> matrices_;
- std::shared_ptr<Assembler> assembler_;
- Matrices<Matrix,1> matrices_;
+ ExternalForce externalForce_;
- VectorType gravityFunctional_;
- std::shared_ptr<const EnergyNorm<ScalarMatrix, ScalarVector>> stateEnergyNorm_;
+ 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
+ // 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_;
- std::vector<std::shared_ptr<GlobalFriction<Matrix, VectorType>>> globalFriction_;
+ // 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);
-
- const std::shared_ptr<Assembler>& assembler() const {
- return assembler_;
- }
-
- const std::shared_ptr<BodyData<dims>>& data() const {
- return bodyData_;
- }
-
- void addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition) {
- leafBoundaryConditions_.push_back(leafBoundaryCondition);
- }
-
- const std::vector<std::shared_ptr<LeafBoundaryCondition>>& leafBoundaryConditions() const {
- return leafBoundaryConditions_;
- }
-
- void leafBoundaryConditions(const std::string& tag, std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const {
- selectedConditions.resize(0);
- for (size_t i=0; i<leafBoundaryConditions_.size(); i++) {
- if (leafBoundaryConditions_[i]->tag() == tag)
- selectedConditions.push_back(leafBoundaryConditions_[i]);
- }
- }
-
- void addBoundaryCondition(std::shared_ptr<FrictionBoundaryCondition> frictionBoundaryCondition) {
- frictionBoundaryConditions_.push_back(frictionBoundaryCondition);
- }
-
- const std::vector<std::shared_ptr<FrictionBoundaryCondition>>& frictionBoundaryConditions() const {
- return frictionBoundaryConditions_;
- }
-
- void addBoundaryCondition(std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition, int level) {
- levelBoundaryConditions_[level].push_back(levelBoundaryCondition);
- }
-
- const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& levelBoundaryConditions() const {
- return levelBoundaryConditions_;
- }
-
- void levelBoundaryConditions(const std::string& tag, std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const {
- selectedConditions.resize(levelBoundaryConditions_.size());
- for (size_t level=0; level<levelBoundaryConditions_.size(); level++) {
- const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
- for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
- if (levelBoundaryConditions[i]->tag() == tag)
- selectedConditions[level].push_back(levelBoundaryConditions[i]);
- }
- }
- }
-
- const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions(int level) const {
- return levelBoundaryConditions_[level];
- }
-
- void levelBoundaryConditions(const std::string& tag, int level, std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const {
- selectedConditions.resize(0);
- const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
- for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
- if (levelBoundaryConditions[i]->tag() == tag)
- selectedConditions.push_back(levelBoundaryConditions[i]);
- }
- }
-
- const DeformedGridType* deformedGrid() const {
- return deformedGrid_.get();
- }
-
- DeformationFunction& deformation() {
- return *deformation_;
- }
-
- void setDeformation(const VectorType& def) {
- deformation_->setDeformation(def);
- }
-
- const DeformedLeafGridView& leafView() const {
- return leafView_;
- }
-
- int leafVertexCount() const {
- return leafVertexCount_;
- }
-
- const Matrices<Matrix,1>& matrices() const {
- return matrices_;
- }
-
- void externalForce(double relativeTime, VectorType& ell) const;
+ 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>>&;
+
+ void setDeformation(const VectorType& def);
+ auto deformation() const -> DeformationFunction&;
+ auto deformedGrid() const -> const DeformedGridType*;
+
+ auto leafView() const -> const DeformedLeafGridView&;
+ auto leafVertexCount() const -> int;
+
+ auto assembler() const -> const std::shared_ptr<Assembler>&;
+ auto matrices() const -> const Matrices<Matrix,1>&;
+
+ 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,
+ std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const;
+ auto leafBoundaryConditions() const -> const std::vector<std::shared_ptr<LeafBoundaryCondition>>&;
+
+ void addBoundaryCondition(std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition, int level);
+ void levelBoundaryConditions(
+ const std::string& tag,
+ std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const;
+ void levelBoundaryConditions(
+ const std::string& tag,
+ int level,
+ std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const;
+ 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;
+ void boundaryNodes(
+ const std::string& tag,
+ BoundaryNodes& nodes) const;
+ void boundaryFunctions(
+ const std::string& tag,
+ BoundaryFunctions& functions) const;
+ void boundaryPatches(
+ const std::string& tag,
+ BoundaryPatches& patches) const;
};
#endif
diff --git a/src/data-structures/contactnetwork.cc b/src/data-structures/contactnetwork.cc
index 93708d26cf42da834ba82861cf69a8570c80ce23..ac118a1ec8f5e443e70bac28a756c84a5491cd03 100644
--- a/src/data-structures/contactnetwork.cc
+++ b/src/data-structures/contactnetwork.cc
@@ -2,124 +2,30 @@
#include "config.h"
#endif
-#include <dune/contact/projections/normalprojection.hh>
#include "contactnetwork.hh"
-template <class GridType, int dimension>
-ContactNetwork<GridType, dimension>::ContactNetwork(int nBodies, int nCouplings, int level) :
- level_(level),
- bodies_(nBodies),
- couplings_(nCouplings),
- deformedGrids_(nBodies),
- leafViews_(nBodies),
- levelViews_(nBodies),
- leafVertexCounts_(nBodies),
- nBodyAssembler_(nBodies, nCouplings)
-{}
-
-template <class GridType, int dimension>
-void ContactNetwork<GridType, dimension>::assemble() {
- for (size_t i=0; i<nBodies(); i++) {
- bodies_[i]->assemble();
- }
-
- totalNodes("dirichlet", totalDirichletNodes_);
-
- // set up dune-contact nBodyAssembler
- nBodyAssembler_.setGrids(deformedGrids_);
-
- for (size_t i=0; i<nCouplings(); i++) {
- nBodyAssembler_.setCoupling(*couplings_[i], i);
- }
-
- nBodyAssembler_.assembleTransferOperator();
- nBodyAssembler_.assembleObstacle();
-}
-
-template <class GridType, int dimension>
-void ContactNetwork<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]);
- }
-}
-
-template <class GridType, int dimension>
-void ContactNetwork<GridType, dimension>::setBodies(const std::vector<std::shared_ptr<Body>>& bodies) {
- assert(nBodies()==bodies.size());
- bodies_ = bodies;
-
- matrices_.elasticity.resize(nBodies());
- matrices_.damping.resize(nBodies());
- matrices_.mass.resize(nBodies());
-
- for (size_t i=0; i<nBodies(); i++) {
- deformedGrids_[i] = bodies_[i]->deformedGrid();
-
- leafViews_[i] = std::make_unique<DeformedLeafGridView>(deformedGrids_[i]->leafGridView());
- levelViews_[i] = std::make_unique<DeformedLevelGridView>(deformedGrids_[i]->levelGridView(0));
-
- leafVertexCounts_[i] = leafViews_[i]->size(dimension);
-
- matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
- matrices_.damping[i] = bodies_[i]->matrices().damping;
- matrices_.mass[i] = bodies_[i]->matrices().mass;
+template <class GridType, int dim>
+void ContactNetwork<GridType, dim>::assemble() {
+ for (size_t i=0; i<levelContactNetworks_.size(); i++) {
+ levelContactNetworks_[i]->assemble();
}
}
-template <class GridType, int dimension>
-void ContactNetwork<GridType, dimension>::setCouplings(const std::vector<std::shared_ptr<CouplingPair>>& couplings) {
- assert(this->nCouplings()==couplings.size());
- couplings_ = couplings;
-}
+template <class GridType, int dim>
+void ContactNetwork<GridType, dim>::addLevelContactNetwork(std::shared_ptr<LevelContactNetwork> levelContactNetwork) {
+ size_t level = levelContactNetwork->level();
+ assert(level < this->size());
-template <class GridType, int dimension>
-void ContactNetwork<GridType, dimension>::constructBody(const std::shared_ptr<BodyData<dimension>>& bodyData,
- const std::shared_ptr<GridType>& grid,
- std::shared_ptr<Body>& body) const
-{
- body = std::make_shared<Body>(bodyData, grid);
+ levelContactNetworks_[level] = levelContactNetwork;
}
-template <class GridType, int dimension>
-void ContactNetwork<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
-{
- coupling = std::make_shared<CouplingPair>();
-
- auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<DeformedLeafBoundaryPatch>>();
- std::shared_ptr<typename CouplingPair::BackEndType> backend = nullptr;
- coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
-}
+template <class GridType, int dim>
+const std::shared_ptr<typename ContactNetwork<GridType, dim>::LevelContactNetwork>& ContactNetwork<GridType, dim>::addLevelContactNetwork(int nBodies, int nCouplings, int level) {
+ assert(level < this->size());
-// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
-template <class GridType, int dimension>
-void ContactNetwork<GridType, dimension>::totalNodes(const std::string& tag, Dune::BitSetVector<dimension>& totalNodes) {
-
- int totalSize = 0;
- for (size_t i=0; i<nBodies(); i++) {
- totalSize += this->body(i)->leafVertexCount();
- }
-
- totalNodes.resize(totalSize);
-
- int idx=0;
- for (size_t i=0; i<nBodies(); i++) {
- const auto& body = this->body(i);
- std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
- body->leafBoundaryConditions(tag, boundaryConditions);
-
- const int idxBackup = idx;
- for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
- const auto& nodes = boundaryConditions[bc]->boundaryNodes();
- for (size_t j=0; j<nodes->size(); j++, idx++)
- for (int k=0; k<dimension; k++)
- totalNodes[idx][k] = (*nodes)[j][k];
- idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
- }
- }
+ levelContactNetworks_[level] = std::make_shared<LevelContactNetwork>(nBodies, nCouplings, level);
+ return levelContactNetworks_[level];
}
#include "contactnetwork_tmpl.cc"
diff --git a/src/data-structures/contactnetwork.hh b/src/data-structures/contactnetwork.hh
index 1a29a90b8eb61004e865194306830f09f1365042..25baf45b74620318079a66894fac23b4703a6948 100644
--- a/src/data-structures/contactnetwork.hh
+++ b/src/data-structures/contactnetwork.hh
@@ -1,145 +1,37 @@
-#ifndef SRC_LEVELCONTACTNETWORK_HH
-#define SRC_LEVELCONTACTNETWORK_HH
+#ifndef SRC_CONTACTNETWORK_HH
+#define SRC_CONTACTNETWORK_HH
-#include <dune/common/fvector.hh>
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/istl/bvector.hh>
-
-#include <dune/contact/common/deformedcontinuacomplex.hh>
-#include <dune/contact/common/couplingpair.hh>
-#include <dune/contact/assemblers/nbodyassembler.hh>
-//#include <dune/contact/assemblers/dualmortarcoupling.hh>
-#include <dune/contact/projections/normalprojection.hh>
-
-#include <dune/tectonic/bodydata.hh>
-#include <dune/tectonic/globalfriction.hh>
-
-#include "../assemblers.hh"
-#include "body.hh"
-#include "enums.hh"
-#include "matrices.hh"
-//#include "multi-body-problem-data/myglobalfrictiondata.hh"
-
-template <class GridType, int dimension> class ContactNetwork {
- private:
- using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+#include "levelcontactnetwork.hh"
+template <class GridType, int dim> class ContactNetwork {
public:
- using Body = Body<GridType, Vector, dimension>;
-
- using DeformedGridType = typename Body::DeformedGridType;
- using DeformedLeafGridView = typename Body::DeformedLeafGridView;
- using DeformedLevelGridView = typename Body::DeformedLevelGridView;
-
- using DeformedLeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
- using DeformedLevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
-
-
-
-
- using Assembler = typename Body::Assembler;
- using Matrix = typename Assembler::Matrix;
- using LocalVector = typename Vector::block_type;
- using ScalarMatrix = typename Assembler::ScalarMatrix;
- using ScalarVector = typename Assembler::ScalarVector;
- using field_type = typename Matrix::field_type;
-
- using CouplingPair = Dune::Contact::CouplingPair<DeformedGridType, DeformedGridType, field_type>;
-
- using Function = Dune::VirtualFunction<double, double>;
+ using LevelContactNetwork = LevelContactNetwork<GridType, dim>;
private:
- const int level_;
-
- std::vector<std::shared_ptr<Body>> bodies_;
- std::vector<std::shared_ptr<CouplingPair>> couplings_;
-
- std::vector<const DeformedGridType*> deformedGrids_;
- std::vector<std::unique_ptr<const DeformedLeafGridView>> leafViews_;
- std::vector<std::unique_ptr<const DeformedLevelGridView>> levelViews_;
- std::vector<size_t> leafVertexCounts_;
-
- Dune::Contact::NBodyAssembler<DeformedGridType, Vector> nBodyAssembler_;
- Dune::BitSetVector<dimension> totalDirichletNodes_;
+ std::vector<std::shared_ptr<LevelContactNetwork>> levelContactNetworks_;
- Matrices<Matrix,2> matrices_;
- // std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionData_;
- std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction_;
+ public:
+ ContactNetwork() {}
- public:
- ContactNetwork(int nBodies, int nCouplings, int level = 0);
+ void addLevelContactNetwork(std::shared_ptr<LevelContactNetwork> levelContactNetwork);
+ const std::shared_ptr<LevelContactNetwork>& addLevelContactNetwork(int nBodies, int nCouplings, int level);
void assemble();
- void assembleFriction(const Config::FrictionModel& frictionModel, const std::vector<ScalarVector>& weightedNormalStress);
-
- int level() const {
- return level_;
- }
-
- void setBodies(const std::vector<std::shared_ptr<Body>>& bodies);
- void setCouplings(const std::vector<std::shared_ptr<CouplingPair>>& couplings);
-
- void constructBody(const std::shared_ptr<BodyData<dimension>>& bodyData,
- const std::shared_ptr<GridType>& grid,
- std::shared_ptr<Body>& body) const;
- void constructCoupling(int nonMortarBodyIdx, int mortarBodyIdx,
- const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
- const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
- std::shared_ptr<CouplingPair>& coupling) const;
-
- void totalNodes(const std::string& tag, Dune::BitSetVector<dimension>& totalNodes);
-
- size_t nBodies() const {
- return bodies_.size();
- }
-
- size_t nCouplings() const {
- return couplings_.size();
- }
-
- const std::vector<size_t>& leafVertexCounts() const {
- return leafVertexCounts_;
- }
- const std::shared_ptr<Body>& body(int i) const {
- return bodies_[i];
+ const std::shared_ptr<LevelContactNetwork>& levelContactNetwork(size_t level) const {
+ return levelContactNetworks_[level];
}
- const std::vector<const DeformedGridType*>& deformedGrids() const {
- return deformedGrids_;
+ void resize(size_t size) {
+ levelContactNetworks_.resize(size);
}
- const DeformedLeafGridView& leafView(int i) const {
- return *leafViews_[i];
+ size_t size() const {
+ return levelContactNetworks_.size();
}
- const DeformedLevelGridView& levelView(int i) const {
- return *levelViews_[i];
- }
-
- int leafVertexCount(int i) const {
- return leafVertexCounts_[i];
- }
-
- Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& nBodyAssembler() {
- return nBodyAssembler_;
- }
-
- const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& nBodyAssembler() const {
- return nBodyAssembler_;
- }
-
- const Matrices<Matrix,2>& matrices() const {
- return matrices_;
- }
-
- const std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>>& globalFriction() const {
- return globalFriction_;
- }
-
- const Dune::BitSetVector<dimension>& totalDirichletNodes() const {
- return totalDirichletNodes_;
- }
+ size_t maxLevel() const {
+ return size()-1;
+ }
};
#endif
diff --git a/src/data-structures/enumparser.hh b/src/data-structures/enumparser.hh
index 34f1c5830f6063d144d6a58dfd092b8cf2896122..47a73f78d0b82e304ef55e99b602d57f3ba0283e 100644
--- a/src/data-structures/enumparser.hh
+++ b/src/data-structures/enumparser.hh
@@ -1,5 +1,5 @@
-#ifndef SRC_ENUMPARSER_HH
-#define SRC_ENUMPARSER_HH
+#ifndef SRC_DATA_STRUCTURES_ENUMPARSER_HH
+#define SRC_DATA_STRUCTURES_ENUMPARSER_HH
// Copyright Carsten Graeser 2012
diff --git a/src/data-structures/globalfrictioncontainer.cc b/src/data-structures/globalfrictioncontainer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..37843c46c76e6c6d0a56c3f6bcba3c6ef9d62af3
--- /dev/null
+++ b/src/data-structures/globalfrictioncontainer.cc
@@ -0,0 +1,114 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+/*
+#include "globalfrictioncontainer.hh"
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i)
+-> IndexObject& {
+ return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i) const
+-> const IndexObject& {
+ return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::size() const
+-> size_t {
+ return globalFriction_.size();
+}
+
+template <class BaseGlobalFriction, size_t depth>
+void GlobalFrictionContainer<BaseGlobalFriction, depth>::resize(std::list<size_t> list) {
+ assert(list.size() <= depth);
+
+ if (list.size() == 0) {
+ globalFriction_.resize(0);
+ } else {
+ globalFriction_.resize(list.front());
+ list.pop_front();
+
+ for (size_t i=0; i<size(); i++) {
+ globalFriction_[i].resize(list);
+ }
+ }
+}
+
+template <class BaseGlobalFriction, size_t depth>
+template <class VectorContainer>
+void GlobalFrictionContainer<BaseGlobalFriction, depth>::updateAlpha(const VectorContainer& newAlpha) {
+ assert(newAlpha.size() == size());
+
+ for (size_t i=0; i<size(); i++) {
+ globalFriction_[i].updateAlpha(newAlpha[i]);
+ }
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction()
+-> GlobalFriction& {
+ return globalFriction_;
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction() const
+-> const GlobalFriction& {
+ return globalFriction_;
+}
+
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i)
+-> IndexObject& {
+ return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i) const
+-> const IndexObject& {
+ return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::size() const
+-> size_t {
+ return globalFriction_.size();
+}
+
+template <class BaseGlobalFriction>
+void GlobalFrictionContainer<BaseGlobalFriction, 1>::resize(std::list<size_t> newSize) {
+ if (newSize.size() > 0) {
+ globalFriction_.resize(newSize.front(), nullptr);
+ } else {
+ globalFriction_.resize(0);
+ }
+}
+
+template <class BaseGlobalFriction>
+template <class Vector>
+void GlobalFrictionContainer<BaseGlobalFriction, 1>::updateAlpha(const Vector& newAlpha) {
+ for (size_t i=0; i<size(); i++) {
+ globalFriction_[i]->updateAlpha(newAlpha);
+ }
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction()
+-> GlobalFriction& {
+ return globalFriction_;
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction() const
+-> const GlobalFriction& {
+ return globalFriction_;
+}
+
+
+#include "globalfrictioncontainer_tmpl.cc"
+*/
diff --git a/src/data-structures/globalfrictioncontainer.hh b/src/data-structures/globalfrictioncontainer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fd2c058b750c041afb01461a313cff9fa5f28979
--- /dev/null
+++ b/src/data-structures/globalfrictioncontainer.hh
@@ -0,0 +1,114 @@
+#ifndef SRC_DATA_STRUCTURES_GLOBALFRICTIONCONTAINER_HH
+#define SRC_DATA_STRUCTURES_GLOBALFRICTIONCONTAINER_HH
+
+#include <vector>
+#include <cstddef>
+#include <list>
+#include <cassert>
+#include <memory>
+
+template <class BaseGlobalFriction, size_t depth>
+class GlobalFrictionContainer {
+public:
+ using IndexObject = GlobalFrictionContainer<BaseGlobalFriction, depth-1>;
+ using GlobalFriction = std::vector<IndexObject>;
+
+ GlobalFrictionContainer() {}
+
+ auto operator[](std::size_t i) -> IndexObject& {
+ return globalFriction_[i];
+ }
+
+ auto operator[](std::size_t i) const -> const IndexObject& {
+ return globalFriction_[i];
+ }
+
+ auto size() const -> size_t {
+ return globalFriction_.size();
+ }
+
+ void resize(std::list<size_t> list) {
+ assert(list.size() <= depth);
+
+ if (list.size() == 0) {
+ globalFriction_.resize(0);
+ } else {
+ globalFriction_.resize(list.front());
+ list.pop_front();
+
+ for (size_t i=0; i<size(); i++) {
+ globalFriction_[i].resize(list);
+ }
+ }
+ }
+
+ template <class VectorContainer>
+ void updateAlpha(const VectorContainer& newAlpha) {
+ assert(newAlpha.size() == size());
+
+ for (size_t i=0; i<size(); i++) {
+ globalFriction_[i].updateAlpha(newAlpha[i]);
+ }
+ }
+
+ auto globalFriction() -> GlobalFriction& {
+ return globalFriction_;
+ }
+
+ auto globalFriction() const -> const GlobalFriction& {
+ return globalFriction_;
+ }
+
+private:
+ GlobalFriction globalFriction_;
+};
+
+
+template <class BaseGlobalFriction>
+class GlobalFrictionContainer<BaseGlobalFriction, 1> {
+public:
+ using IndexObject = std::shared_ptr<BaseGlobalFriction>;
+ using GlobalFriction = std::vector<IndexObject>;
+
+ GlobalFrictionContainer() {}
+
+ auto operator[](std::size_t i) -> IndexObject& {
+ return globalFriction_[i];
+ }
+
+ auto operator[](std::size_t i) const -> const IndexObject& {
+ return globalFriction_[i];
+ }
+
+ auto size() const -> size_t {
+ return globalFriction_.size();
+ }
+
+ void resize(std::list<size_t> newSize) {
+ if (newSize.size() > 0) {
+ globalFriction_.resize(newSize.front(), nullptr);
+ } else {
+ globalFriction_.resize(0);
+ }
+ }
+
+ template <class Vector>
+ void updateAlpha(const Vector& newAlpha) {
+ for (size_t i=0; i<size(); i++) {
+ globalFriction_[i]->updateAlpha(newAlpha);
+ }
+ }
+
+ auto globalFriction() -> GlobalFriction& {
+ return globalFriction_;
+ }
+
+ auto globalFriction() const -> const GlobalFriction& {
+ return globalFriction_;
+ }
+
+private:
+ GlobalFriction globalFriction_;
+};
+
+#endif
diff --git a/src/data-structures/levelcontactnetwork.cc b/src/data-structures/levelcontactnetwork.cc
index 8f58d4d250fa29abe3c268667db81a1fd9da675a..bf3ca67689f6f15e418261f185462aef33d661ab 100644
--- a/src/data-structures/levelcontactnetwork.cc
+++ b/src/data-structures/levelcontactnetwork.cc
@@ -6,8 +6,14 @@
#include "levelcontactnetwork.hh"
+#include "../utils/tobool.hh"
+#include "../utils/debugutils.hh"
+
template <class GridType, int dimension>
-LevelContactNetwork<GridType, dimension>::LevelContactNetwork(int nBodies, int nCouplings, int level) :
+LevelContactNetwork<GridType, dimension>::LevelContactNetwork(
+ int nBodies,
+ int nCouplings,
+ int level) :
level_(level),
bodies_(nBodies),
couplings_(nCouplings),
@@ -22,10 +28,9 @@ template <class GridType, int dimension>
void LevelContactNetwork<GridType, dimension>::assemble() {
for (size_t i=0; i<nBodies(); i++) {
bodies_[i]->assemble();
+ //printBasisDofLocation(bodies_[i]->assembler()->vertexBasis); //TODO remove after debugging
}
- totalNodes("dirichlet", totalDirichletNodes_);
-
// set up dune-contact nBodyAssembler
nBodyAssembler_.setGrids(deformedGrids_);
@@ -38,7 +43,10 @@ void LevelContactNetwork<GridType, dimension>::assemble() {
}
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::assembleFriction(const Config::FrictionModel& frictionModel, const std::vector<ScalarVector>& weightedNormalStress) {
+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]);
}
@@ -57,7 +65,11 @@ void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::
deformedGrids_[i] = bodies_[i]->deformedGrid();
leafViews_[i] = std::make_unique<DeformedLeafGridView>(deformedGrids_[i]->leafGridView());
- levelViews_[i] = std::make_unique<DeformedLevelGridView>(deformedGrids_[i]->levelGridView(0));
+
+ levelViews_[i].resize(deformedGrids_[i]->maxLevel()+1);
+ for (size_t j=0; j<levelViews_[i].size(); j++) {
+ levelViews_[i][j] = std::make_unique<DeformedLevelGridView>(deformedGrids_[i]->levelGridView(j));
+ }
leafVertexCounts_[i] = leafViews_[i]->size(dimension);
@@ -74,19 +86,21 @@ void LevelContactNetwork<GridType, dimension>::setCouplings(const std::vector<st
}
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::constructBody(const std::shared_ptr<BodyData<dimension>>& bodyData,
- const std::shared_ptr<GridType>& grid,
- std::shared_ptr<Body>& body) const
-{
+void LevelContactNetwork<GridType, dimension>::constructBody(
+ const std::shared_ptr<BodyData<dimension>>& bodyData,
+ const std::shared_ptr<GridType>& grid,
+ std::shared_ptr<Body>& body) const {
+
body = std::make_shared<Body>(bodyData, grid);
}
template <class GridType, int dimension>
-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
-{
+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 {
+
coupling = std::make_shared<CouplingPair>();
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<DeformedLeafBoundaryPatch>>();
@@ -96,7 +110,9 @@ void LevelContactNetwork<GridType, dimension>::constructCoupling(int nonMortarBo
// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::totalNodes(const std::string& tag, Dune::BitSetVector<dimension>& totalNodes) {
+void LevelContactNetwork<GridType, dimension>::totalNodes(
+ const std::string& tag,
+ Dune::BitSetVector<dimension>& totalNodes) const {
int totalSize = 0;
for (size_t i=0; i<nBodies(); i++) {
@@ -122,4 +138,174 @@ void LevelContactNetwork<GridType, dimension>::totalNodes(const std::string& tag
}
}
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::boundaryPatchNodes(
+ const std::string& tag,
+ BoundaryPatchNodes& nodes) const {
+
+ nodes.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ this->body(i)->boundaryPatchNodes(tag, nodes[i]);
+ }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::boundaryNodes(
+ const std::string& tag,
+ BoundaryNodes& nodes) const {
+
+ nodes.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ this->body(i)->boundaryNodes(tag, nodes[i]);
+ }
+}
+
+// collects all leaf boundary functions from the different bodies identified by "tag"
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::boundaryFunctions(
+ const std::string& tag,
+ BoundaryFunctions& functions) const {
+
+ functions.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ this->body(i)->boundaryFunctions(tag, functions[i]);
+ }
+}
+
+// collects all leaf boundary patches from the different bodies identified by "tag"
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::boundaryPatches(
+ const std::string& tag,
+ BoundaryPatches& patches) const {
+
+ patches.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ this->body(i)->boundaryPatches(tag, patches[i]);
+ }
+}
+
+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();
+ }
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::level() const
+-> int {
+
+ return level_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nBodies() const
+-> size_t {
+
+ return bodies_.size();
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nCouplings() const
+-> size_t {
+
+ return couplings_.size();
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::leafVertexCounts() const
+-> const std::vector<size_t>& {
+
+ return leafVertexCounts_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::body(int i) const
+-> const std::shared_ptr<Body>& {
+
+ return bodies_[i];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::deformedGrids() const
+-> const std::vector<const DeformedGridType*>& {
+
+ return deformedGrids_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::leafView(int bodyID) const
+-> const DeformedLeafGridView& {
+
+ return *leafViews_[bodyID];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::levelViews(int bodyID) const
+-> const std::vector<std::unique_ptr<const DeformedLevelGridView>>& {
+
+ return levelViews_[bodyID];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::levelView(int bodyID, int level) const
+-> const DeformedLevelGridView& {
+
+ return *levelViews_[bodyID][level];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::leafVertexCount(int bodyID) const
+-> int {
+
+ return leafVertexCounts_[bodyID];
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nBodyAssembler()
+-> Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& {
+
+ return nBodyAssembler_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::nBodyAssembler() const
+-> const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& {
+
+ return nBodyAssembler_;
+}
+
+template <class GridType, int dimension>
+auto LevelContactNetwork<GridType, dimension>::matrices() const
+-> const Matrices<Matrix,2>& {
+
+ return matrices_;
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::globalFriction(GlobalFrictionContainer& globalFriction) const {
+ globalFriction.resize({nBodies()});
+
+ for (size_t i=0; i<nBodies(); i++) {
+ globalFriction[i] = bodies_[i]->globalFriction();
+ }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::externalForces(ExternalForces& externalForces) const {
+ externalForces.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ externalForces[i] = &bodies_[i]->externalForce();
+ }
+}
+
+
#include "levelcontactnetwork_tmpl.cc"
diff --git a/src/data-structures/levelcontactnetwork.hh b/src/data-structures/levelcontactnetwork.hh
index ea47b6efdd690c45732ab4096ec4c09483f10c13..f6c1bf194f3889ff2effa0a446f317df4e09621f 100644
--- a/src/data-structures/levelcontactnetwork.hh
+++ b/src/data-structures/levelcontactnetwork.hh
@@ -1,5 +1,5 @@
-#ifndef SRC_LEVELCONTACTNETWORK_HH
-#define SRC_LEVELCONTACTNETWORK_HH
+#ifndef SRC_DATA_STRUCTURES_LEVELCONTACTNETWORK_HH
+#define SRC_DATA_STRUCTURES_LEVELCONTACTNETWORK_HH
#include <dune/common/fvector.hh>
#include <dune/common/bitsetvector.hh>
@@ -13,133 +13,128 @@
#include <dune/contact/projections/normalprojection.hh>
#include <dune/tectonic/bodydata.hh>
-#include <dune/tectonic/globalfriction.hh>
#include "../assemblers.hh"
#include "body.hh"
#include "enums.hh"
+#include "globalfrictioncontainer.hh"
#include "matrices.hh"
//#include "multi-body-problem-data/myglobalfrictiondata.hh"
template <class GridType, int dimension> class LevelContactNetwork {
- private:
- using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+ private:
+ using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
- public:
- using Body = Body<GridType, Vector, dimension>;
+ public:
+ using Body = Body<GridType, Vector, dimension>;
- using DeformedGridType = typename Body::DeformedGridType;
- using DeformedLeafGridView = typename Body::DeformedLeafGridView;
- using DeformedLevelGridView = typename Body::DeformedLevelGridView;
+ using DeformedGridType = typename Body::DeformedGridType;
+ using DeformedLeafGridView = typename Body::DeformedLeafGridView;
+ using DeformedLevelGridView = typename Body::DeformedLevelGridView;
- using DeformedLeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
- using DeformedLevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
+ using DeformedLeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
+ using DeformedLevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
+ using Assembler = typename Body::Assembler;
+ using Matrix = typename Assembler::Matrix;
+ using LocalVector = typename Vector::block_type;
+ using ScalarMatrix = typename Assembler::ScalarMatrix;
+ using ScalarVector = typename Assembler::ScalarVector;
+ using field_type = typename Matrix::field_type;
+
+ using CouplingPair = Dune::Contact::CouplingPair<DeformedGridType, DeformedGridType, field_type>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+ using BoundaryFunctions = std::vector<typename Body::BoundaryFunctions>;
+ using BoundaryNodes = std::vector<typename Body::BoundaryNodes>;
+ using BoundaryPatchNodes = std::vector<typename Body::BoundaryPatchNodes>;
+ using BoundaryPatches = std::vector<typename Body::BoundaryPatches>;
+ using GlobalFrictionContainer = GlobalFrictionContainer<typename Body::GlobalFriction, 2>;
- using Assembler = typename Body::Assembler;
- using Matrix = typename Assembler::Matrix;
- using LocalVector = typename Vector::block_type;
- using ScalarMatrix = typename Assembler::ScalarMatrix;
- using ScalarVector = typename Assembler::ScalarVector;
- using field_type = typename Matrix::field_type;
-
- using CouplingPair = Dune::Contact::CouplingPair<DeformedGridType, DeformedGridType, field_type>;
-
- using Function = Dune::VirtualFunction<double, double>;
-
- private:
- const int level_;
-
- std::vector<std::shared_ptr<Body>> bodies_;
- std::vector<std::shared_ptr<CouplingPair>> couplings_;
-
- std::vector<const DeformedGridType*> deformedGrids_;
- std::vector<std::unique_ptr<const DeformedLeafGridView>> leafViews_;
- std::vector<std::unique_ptr<const DeformedLevelGridView>> levelViews_;
- std::vector<size_t> leafVertexCounts_;
-
- Dune::Contact::NBodyAssembler<DeformedGridType, Vector> nBodyAssembler_;
- Dune::BitSetVector<dimension> totalDirichletNodes_;
-
- Matrices<Matrix,2> matrices_;
- // std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionData_;
- std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction_;
-
- public:
+ using StateEnergyNorms = std::vector<const typename Body::StateEnergyNorm *>;
+ using ExternalForces = std::vector<const typename Body::ExternalForce *>;
+
+ public:
LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
void assemble();
- void assembleFriction(const Config::FrictionModel& frictionModel, const std::vector<ScalarVector>& weightedNormalStress);
-
- int level() const {
- return level_;
- }
+ 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 constructBody(const std::shared_ptr<BodyData<dimension>>& bodyData,
- const std::shared_ptr<GridType>& grid,
- std::shared_ptr<Body>& body) const;
- void constructCoupling(int nonMortarBodyIdx, int mortarBodyIdx,
- const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
- const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
- std::shared_ptr<CouplingPair>& coupling) const;
+ void constructBody(
+ const std::shared_ptr<BodyData<dimension>>& bodyData,
+ const std::shared_ptr<GridType>& grid,
+ std::shared_ptr<Body>& body) const;
+ void constructCoupling(
+ int nonMortarBodyIdx,
+ int mortarBodyIdx,
+ const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
+ const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
+ std::shared_ptr<CouplingPair>& coupling) const;
+
+ // getter
+ void totalNodes(
+ const std::string& tag,
+ Dune::BitSetVector<dimension>& totalNodes) const;
+ void boundaryPatches(
+ const std::string& tag,
+ BoundaryPatches& patches) const;
+ void boundaryPatchNodes(
+ const std::string& tag,
+ BoundaryPatchNodes& nodes) const;
+ void boundaryNodes(
+ const std::string& tag,
+ BoundaryNodes& nodes) const;
+ void boundaryFunctions(
+ const std::string& tag,
+ BoundaryFunctions& functions) const;
+
+ void stateEnergyNorms(StateEnergyNorms& stateEnergyNorms) const;
+ void globalFriction(GlobalFrictionContainer& globalFriction) const;
+ void externalForces(ExternalForces& externalForces) const;
+
+ auto level() const -> int;
- void totalNodes(const std::string& tag, Dune::BitSetVector<dimension>& totalNodes);
+ auto nBodies() const -> size_t;
+ auto nCouplings() const -> size_t;
- size_t nBodies() const {
- return bodies_.size();
- }
+ auto body(int i) const -> const std::shared_ptr<Body>&;
- size_t nCouplings() const {
- return couplings_.size();
- }
+ auto deformedGrids() const -> const std::vector<const DeformedGridType*>& ;
- const std::vector<size_t>& leafVertexCounts() const {
- return leafVertexCounts_;
- }
+ auto leafVertexCounts() const -> const std::vector<size_t>&;
+ auto leafVertexCount(int bodyID) const -> int;
- const std::shared_ptr<Body>& body(int i) const {
- return bodies_[i];
- }
+ auto leafView(int bodyID) const -> const DeformedLeafGridView&;
+ auto levelView(int bodyID, int level = 0) const -> const DeformedLevelGridView&;
+ auto levelViews(int bodyID) const -> const std::vector<std::unique_ptr<const DeformedLevelGridView>>&;
- const std::vector<const DeformedGridType*>& deformedGrids() const {
- return deformedGrids_;
- }
+ auto nBodyAssembler() -> Dune::Contact::NBodyAssembler<DeformedGridType, Vector>&;
+ auto nBodyAssembler() const -> const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>&;
- const DeformedLeafGridView& leafView(int i) const {
- return *leafViews_[i];
- }
+ auto matrices() const -> const Matrices<Matrix,2>&;
- const DeformedLevelGridView& levelView(int i) const {
- return *levelViews_[i];
- }
- int leafVertexCount(int i) const {
- return leafVertexCounts_[i];
- }
- Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& nBodyAssembler() {
- return nBodyAssembler_;
- }
+ private:
+ const int level_;
- const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& nBodyAssembler() const {
- return nBodyAssembler_;
- }
+ std::vector<std::shared_ptr<Body>> bodies_;
+ std::vector<std::shared_ptr<CouplingPair>> couplings_;
- const Matrices<Matrix,2>& matrices() const {
- return matrices_;
- }
+ 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_;
- const std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>>& globalFriction() const {
- return globalFriction_;
- }
+ Dune::Contact::NBodyAssembler<DeformedGridType, Vector> nBodyAssembler_;
- const Dune::BitSetVector<dimension>& totalDirichletNodes() const {
- return totalDirichletNodes_;
- }
+ Matrices<Matrix,2> matrices_;
};
#endif
diff --git a/src/data-structures/levelcontactnetwork_tmpl.cc b/src/data-structures/levelcontactnetwork_tmpl.cc
index b0c5200f198c29fe30a2407919ad55b94b36f461..1a4be09fc127345d24d54b09e9baef4f3b1f7e83 100644
--- a/src/data-structures/levelcontactnetwork_tmpl.cc
+++ b/src/data-structures/levelcontactnetwork_tmpl.cc
@@ -3,5 +3,6 @@
#endif
#include "../explicitgrid.hh"
+#include "levelcontactnetwork.hh"
template class LevelContactNetwork<Grid, MY_DIM>;
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 44b84f91086a1f35f01dd2e973c4c11cc3c088ea..275697da77effe110850f77c239a4d1981b9a41c 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -99,7 +99,6 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
std::vector<const Matrix*> matrices_ptr(_matrices.size());
for (size_t i=0; i<matrices_ptr.size(); i++) {
matrices_ptr[i] = _matrices[i].get();
- print(*matrices_ptr[i], "matrix " + std::to_string(i));
}
/*std::vector<Matrix> matrices(velocityMatrices.size());
@@ -129,10 +128,15 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
using LinearFactory = SolverFactory<DeformedGridType, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
LinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
- print(_dirichletNodes, "dirichletNodes: ");
- print(bilinearForm, "matrix: ");
+ /* std::vector<Dune::BitSetVector<dims>> bodyDirichletNodes;
+ nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
+ for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+ print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+ }*/
+
+ /* print(bilinearForm, "matrix: ");
print(totalX, "totalX: ");
- print(totalRhs, "totalRhs: ");
+ print(totalRhs, "totalRhs: ");*/
auto multigridStep = factory.getStep();
multigridStep->setProblem(bilinearForm, totalX, totalRhs);
@@ -163,13 +167,14 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
ell0[i].resize(u[i].size());
ell0[i] = 0.0;
- levelContactNetwork.body(i)->externalForce(relativeTime, ell0[i]);
+ levelContactNetwork.body(i)->externalForce()(relativeTime, ell0[i]);
}
// Initial displacement: Start from a situation of minimal stress,
// which is automatically attained in the case [v = 0 = a].
// Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
- Dune::BitSetVector<dims> dirichletNodes = levelContactNetwork.totalDirichletNodes();
+ Dune::BitSetVector<dims> dirichletNodes;
+ levelContactNetwork.totalNodes("dirichlet", dirichletNodes);
for (size_t i=0; i<dirichletNodes.size(); i++) {
bool val = false;
for (size_t d=0; d<dims; d++) {
@@ -210,7 +215,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
}
- Dune::BitSetVector<dims> noNodes(levelContactNetwork.totalDirichletNodes().size());
+ Dune::BitSetVector<dims> noNodes(dirichletNodes.size(), false);
solveLinearProblem(noNodes, levelContactNetwork.matrices().mass, accelerationRHS, a,
parset.sub("a0.solver"));
}
diff --git a/src/factories/cantorfactory.cc b/src/factories/cantorfactory.cc
index 4c30f26697f55681df4663adafd5238e90e76492..4df27f770be9ceb4e6de4a7454ab847978d2e2f1 100644
--- a/src/factories/cantorfactory.cc
+++ b/src/factories/cantorfactory.cc
@@ -8,69 +8,32 @@
#include <dune/contact/projections/normalprojection.hh>
+#include "../utils/almostequal.hh"
+
#include "cantorfactory.hh"
#include "../multi-body-problem-data/bc.hh"
#include "../multi-body-problem-data/cuboidgeometry.hh"
#include "../multi-body-problem-data/myglobalfrictiondata.hh"
#include "../multi-body-problem-data/weakpatch.hh"
+#include "../multi-body-problem-data/cubegridconstructor.hh"
#include "../utils/diameter.hh"
-template <int dim>
-class Cube {
-private:
- using Vector = Dune::FieldVector<double, dim>;
-
- Vector A_; // two corners that are diagonal define cube in any space dimension
- Vector B_;
-
- bool invariant_; // flag to set if Cube can be split()
-
- /* bool isDiagonal(const Vector& A, const Vector& B) const {
- Vector diff = A;
- diff -= B;
+template <class GridType, int dim>
+void CantorFactory<GridType, dim>::setBodies() {
- for (size_t d=0; d<dim; d++) {
- if (std::abs(diff[d])< 1e-14)
- return false;
- }
- return true;
- } */
-public:
- Cube(bool invariant = false) {
- invariant_ = invariant;
- }
- Cube(const Vector& A, const Vector& B, bool invariant = false) {
- setCorners(A, B);
- invariant_ = invariant;
- }
+ for (size_t level=0; level<=maxLevel_; level++) {
+ const LevelCubes& levelCubes = cubes_[level];
- void setCorners(const Vector& A, const Vector& B) {
- //assert(isDiagonal(corners[0], corners[1]));
- A_ = A;
- B_ = B;
- }
+ for (size_t i=0; i<levelCubes.size(); i++) {
+ const std::shared_ptr<Cube>& cube = levelCubes[i];
- void split(std::vector<Cube>& newCubes) const {
- assert(!invariant_);
-
- Vector midPoint = A_;
- midPoint += 1/2*(B_ - A_);
-
- newCubes.push_back();
+ CubeGridConstructor
+ }
}
-};
-
-template <class GridType, int dims>
-void CantorFactory<GridType, dims>::setBodies() {
- // set up cuboid geometries
-
- double const length = 1.00;
- double const width = 0.27;
- double const weakLength = 0.20;
std::vector<LocalVector> origins(this->bodyCount_);
std::vector<LocalVector> lowerWeakOrigins(this->bodyCount_);
@@ -133,18 +96,6 @@ void CantorFactory<GridType, dims>::setBodies() {
getWeakPatch<LocalVector>(this->parset_.sub("boundary.friction.weakening"), cuboidGeometry, *lowerWeakPatch, *upperWeakPatch);
refine(*grids[i], *lowerWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
refine(*grids[i], *upperWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
-
- // determine minDiameter and maxDiameter
- double minDiameter = std::numeric_limits<double>::infinity();
- double maxDiameter = 0.0;
- for (auto &&e : elements(grids[i]->leafGridView())) {
- auto const geometry = e.geometry();
- auto const diam = diameter(geometry);
- minDiameter = std::min(minDiameter, diam);
- maxDiameter = std::max(maxDiameter, diam);
- }
- std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
- std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
}
for (size_t i=0; i<this->bodyCount_; i++) {
@@ -152,8 +103,8 @@ void CantorFactory<GridType, dims>::setBodies() {
}
}
-template <class GridType, int dims>
-void CantorFactory<GridType, dims>::setCouplings() {
+template <class GridType, int dim>
+void CantorFactory<GridType, dim>::setCouplings() {
const auto deformedGrids = this->contactNetwork_.deformedGrids();
for (size_t i=0; i<this->bodyCount_; i++) {
@@ -175,10 +126,10 @@ void CantorFactory<GridType, dims>::setCouplings() {
}
}
-template <class GridType, int dims>
-void CantorFactory<GridType, dims>::setBoundaryConditions() {
- using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
- using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
+template <class GridType, int dim>
+void CantorFactory<GridType, dim>::setBoundaryConditions() {
+ using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dim>;
+ using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dim>;
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
@@ -214,7 +165,7 @@ void CantorFactory<GridType, dims>::setBoundaryConditions() {
// Dirichlet Boundary
// identify Dirichlet nodes on leaf level
- std::shared_ptr<Dune::BitSetVector<dims>> dirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
+ std::shared_ptr<Dune::BitSetVector<dim>> dirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
for (int j=0; j<leafVertexCount; j++) {
if (leafFaces_[i]->right.containsVertex(j))
(*dirichletNodes)[j][0] = true;
@@ -235,4 +186,219 @@ void CantorFactory<GridType, dims>::setBoundaryConditions() {
}
}
+
+class CantorIndexHierarchy {
+ public:
+ typedef std::map<std::pair<int,int>, bool> LevelCantorIndices;
+
+ private:
+ const size_t maxLevel_;
+ std::vector<size_t> levelN_;
+ std::vector<LevelCantorIndices> cantorIndexHierarchy_;
+
+ void prolongIndices(size_t currentLevel, size_t newLevel) {
+ const LevelCantorIndices& indices = cantorIndexHierarchy_[currentLevel];
+ LevelCantorIndices& newIndices = cantorIndexHierarchy_[newLevel];
+
+ size_t offset = levelN_[currentLevel];
+
+ std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
+ std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
+ for (; it!=endIt; ++it) {
+ int xID = it->first.first;
+ int yID = it->first.second;
+ newIndices[std::make_pair(xID, yID)] = true;
+ newIndices[std::make_pair(xID+offset, yID)] = true;
+ newIndices[std::make_pair(xID, yID+offset)] = true;
+ }
+
+ size_t doubleOffset = 2*offset;
+ for (size_t i=offset+1; i<=doubleOffset; i=i+2) {
+ newIndices[std::make_pair(doubleOffset, i)] = true;
+ newIndices[std::make_pair(i, doubleOffset)] = true;
+ }
+ }
+
+ void print(const LevelCantorIndices& indices) const {
+ std::cout << "LevelCantorIndices: " << std::endl;
+
+ std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
+ std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
+ for (; it!=endIt; ++it) {
+ std::cout << "(" << it->first.first << ", " << it->first.second << ")"<< std::endl;
+ }
+ }
+
+ public:
+ CantorIndexHierarchy(size_t maxLevel) :
+ maxLevel_(maxLevel)
+ {
+ levelN_.resize(maxLevel_+1);
+ cantorIndexHierarchy_.resize(maxLevel_+1);
+
+ // init levelCantorIndices on level 0
+ LevelCantorIndices& initCantorIndices = cantorIndexHierarchy_[0];
+ initCantorIndices[std::make_pair(0, 1)] = true;
+ initCantorIndices[std::make_pair(1, 0)] = true;
+ initCantorIndices[std::make_pair(1, 2)] = true;
+ initCantorIndices[std::make_pair(2, 1)] = true;
+
+ for (size_t i=0; i<maxLevel_; i++) {
+ levelN_[i] = std::pow(2, i+1);
+ prolongIndices(i, i+1);
+ }
+
+ levelN_[maxLevel_] = std::pow(2, maxLevel_+1);
+ }
+
+ const LevelCantorIndices& levelCantorIndices(size_t i) const {
+ return cantorIndexHierarchy_[i];
+ }
+
+ size_t levelN(const size_t i) const {
+ return levelN_[i];
+ }
+};
+
+
+template <class GridType>
+class CantorFaultFactory
+{
+ //! Parameter for mapper class
+ template<int dim>
+ struct FaceMapperLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim-1;
+ }
+ };
+
+ protected:
+ typedef OscUnitCube<GridType, 2> GridOb;
+ static const int dimworld = GridType::dimensionworld;
+ static const int dim = GridType::dimension;
+ typedef typename GridType::ctype ctype;
+ typedef typename GridType::LevelGridView GV;
+
+ typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
+
+ private:
+ const std::map<size_t, size_t>& levelToCantorLevel_;
+ const int coarseResolution_;
+ const size_t maxLevel_;
+ const size_t maxCantorLevel_;
+ const CantorIndexHierarchy cantorIndexHierarchy_;
+ const int coarseGridN_;
+
+ std::shared_ptr<GridType> grid_;
+ std::vector<double> levelResolutions_;
+
+ std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;
+
+private:
+
+
+ bool computeID(Dune::FieldVector<ctype, dimworld> vertex, const size_t gridN, std::pair<size_t, size_t>& IDs) const {
+ ctype xID = vertex[0]*gridN;
+ ctype yID = vertex[1]*gridN;
+
+ ctype x = 0;
+ ctype y = 0;
+
+ bool xIsInt = almost_equal(std::modf(xID, &x), 0.0, 2);
+ bool yIsInt = almost_equal(std::modf(yID, &y), 0.0, 2);
+
+ if (xIsInt && yIsInt) {
+ IDs = std::make_pair((size_t) x, (size_t) y);
+ return true;
+ }
+
+ if (xIsInt) {
+ y = std::ceil(yID);
+ if ((size_t) y % 2==0)
+ y = y-1;
+ IDs = std::make_pair((size_t) x, (size_t) y);
+ return true;
+ }
+
+ if (yIsInt) {
+ x = std::ceil(xID);
+ if ((size_t) x % 2==0)
+ x = x-1;
+ IDs = std::make_pair((size_t) x, (size_t) y);
+ return true;
+ }
+
+ return false;
+ }
+
+public:
+ //setup
+ CantorFaultFactory(const std::map<size_t, size_t>& levelToCantorLevel, const int coarseResolution, const size_t maxLevel, const size_t maxCantorLevel) :
+ levelToCantorLevel_(levelToCantorLevel),
+ coarseResolution_(coarseResolution),
+ maxLevel_(maxLevel),
+ maxCantorLevel_(maxCantorLevel),
+ cantorIndexHierarchy_(maxCantorLevel_),
+ coarseGridN_(std::pow(2, coarseResolution_)),
+ interfaceNetwork_(nullptr)
+ {
+ Dune::UGGrid<dim>::setDefaultHeapSize(4000);
+ GridOb unitCube(coarseGridN_);
+ grid_ = unitCube.grid();
+
+ grid_->globalRefine(maxLevel_);
+
+ levelResolutions_.resize(maxLevel_+1);
+
+ // init interface network
+ interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
+
+ for (size_t i=0; i<=maxLevel_; i++) {
+ if (i>0) {
+ interfaceNetwork_->prolongLevel(i-1, i);
+ }
+
+ if (levelToCantorLevel_.count(i)) {
+
+ levelResolutions_[i] = std::pow(2, coarseResolution_+i);
+ const size_t levelResolution = levelResolutions_[i];
+ const size_t cantorResolution = cantorIndexHierarchy_.levelN(levelToCantorLevel_.at(i));
+
+ if (2*levelResolution<cantorResolution)
+ DUNE_THROW(Dune::Exception, "Level " + std::to_string(i) + " does not support cantor set with resolution " + std::to_string(cantorResolution));
+
+ const typename CantorIndexHierarchy::LevelCantorIndices& levelCantorIndices = cantorIndexHierarchy_.levelCantorIndices(levelToCantorLevel_.at(i));
+ std::shared_ptr<LevelInterfaceNetwork<GV>> levelInterfaceNetwork = interfaceNetwork_->levelInterfaceNetworkPtr(i);
+
+ const GV& gridView = levelInterfaceNetwork->levelGridView();
+ FaceMapper intersectionMapper(gridView);
+ std::vector<bool> intersectionHandled(intersectionMapper.size(),false);
+
+ for (const auto& elem:elements(gridView)) {
+ for (const auto& isect:intersections(gridView, elem)) {
+
+ if (intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)])
+ continue;
+
+ intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)]=true;
+
+ if (isect.boundary())
+ continue;
+
+ std::pair<size_t, size_t> intersectionID;
+ bool isAdmissibleID = computeID(isect.geometry().center(), cantorResolution, intersectionID);
+
+ if (isAdmissibleID && levelCantorIndices.count(intersectionID)) {
+ levelInterfaceNetwork->addIntersection(isect, i);
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+
#include "cantorfactory_tmpl.cc"
diff --git a/src/factories/cantorfactory.hh b/src/factories/cantorfactory.hh
index 8818fed94c60ab15da62d023a4e8ee166b404886..7e58c9ad159cdab9168a53eaf62311c55d165578 100644
--- a/src/factories/cantorfactory.hh
+++ b/src/factories/cantorfactory.hh
@@ -10,11 +10,16 @@
#include "contactnetworkfactory.hh"
#include "../multi-body-problem-data/mybody.hh"
-#include "../multi-body-problem-data/mygrids.hh"
+#include "../multi-body-problem-data/grid/mygrids.hh"
-template <class GridType, int dims> class CantorFactory : public ContactNetworkFactory<GridType, dims>{
+
+
+
+template <class GridType, int dim> class CantorFactory : public ContactNetworkFactory<GridType, dim>{
private:
- using Base = ContactNetworkFactory<GridType, dims>;
+ using Base = ContactNetworkFactory<GridType, dim>;
+
+ using LevelCubes = std::vector<std::shared_ptr<Cube>>;
using LocalVector = typename Base::ContactNetwork::LocalVector;
@@ -28,11 +33,11 @@ template <class GridType, int dims> class CantorFactory : public ContactNetworkF
using LevelFaces = MyFaces<DeformedLevelGridView>;
private:
- const std::shared_ptr<MyBodyData<dims>> bodyData_; // material properties of bodies
+ const std::shared_ptr<MyBodyData<dim>> bodyData_; // material properties of bodies
GridsConstructor<GridType>* gridConstructor_;
- std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
+ std::vector<LevelCubes> cubes_;
std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
@@ -46,7 +51,7 @@ template <class GridType, int dims> class CantorFactory : public ContactNetworkF
public:
CantorFactory(const Dune::ParameterTree& parset) :
Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
- bodyData_(std::make_shared<MyBodyData<dims>>(this->parset_, zenith_())),
+ bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_, zenith_())),
cuboidGeometries_(this->bodyCount_),
leafFaces_(this->bodyCount_),
levelFaces_(this->bodyCount_),
diff --git a/src/factories/contactnetworkfactory.hh b/src/factories/contactnetworkfactory.hh
index bee811e7ed717428a93d5c0bdb56cc8bd782401a..90298ad6803b98a8a7cd6053d4b805e3128207bc 100644
--- a/src/factories/contactnetworkfactory.hh
+++ b/src/factories/contactnetworkfactory.hh
@@ -4,22 +4,21 @@
#include <dune/common/parametertree.hh>
#include "../data-structures/contactnetwork.hh"
+#include "../data-structures/levelcontactnetwork.hh"
-template <class GridType, int dims>
+#include "levelcontactnetworkfactory.hh"
+
+template <class GridType, int dim>
class ContactNetworkFactory {
public:
- using ContactNetwork = ContactNetwork<GridType, dims>;
+ using ContactNetwork = ContactNetwork<GridType, dim>;
protected:
- using Body = typename ContactNetwork::Body;
- using CouplingPair = typename ContactNetwork::CouplingPair;
+ using LevelContactNetworkFactory = LevelContactNetworkFactory<GridType, dim>;
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<LevelContactNetworkFactory>> levelContactNetworkFactories_;
ContactNetwork contactNetwork_;
@@ -29,20 +28,16 @@ class ContactNetworkFactory {
virtual void setBoundaryConditions() = 0;
public:
- ContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
+ ContactNetworkFactory(const Dune::ParameterTree& parset) :
parset_(parset),
- bodyCount_(bodyCount),
- couplingCount_(couplingCount),
- bodies_(bodyCount_),
- couplings_(couplingCount_),
- contactNetwork_(bodyCount_, couplingCount_) {}
+ contactNetwork_() {}
void build() {
setBodies();
- contactNetwork_.setBodies(bodies_);
+ contactNetwork_.setBodies();
setCouplings();
- contactNetwork_.setCouplings(couplings_);
+ contactNetwork_.setCouplings();
setBoundaryConditions();
}
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index 504ddff81aa9df3cce6434af4c1e9445890976c5..0891faeed91b05c974393734c48b20442ca01a2c 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -7,7 +7,7 @@
#include <dune/contact/projections/normalprojection.hh>
#include "../multi-body-problem-data/bc.hh"
-#include "../multi-body-problem-data/cuboidgeometry.hh"
+#include "../multi-body-problem-data/grid/cuboidgeometry.hh"
#include "../multi-body-problem-data/myglobalfrictiondata.hh"
#include "../multi-body-problem-data/weakpatch.hh"
diff --git a/src/factories/stackedblocksfactory.hh b/src/factories/stackedblocksfactory.hh
index cf486808766bd42e7b9fdc0c606fa753630ba6f4..ed38c0ff4a01b9564a9e297f4e393d52eb36e83b 100644
--- a/src/factories/stackedblocksfactory.hh
+++ b/src/factories/stackedblocksfactory.hh
@@ -10,7 +10,7 @@
#include "levelcontactnetworkfactory.hh"
#include "../multi-body-problem-data/mybody.hh"
-#include "../multi-body-problem-data/mygrids.hh"
+#include "../multi-body-problem-data/grid/mygrids.hh"
template <class GridType, int dims> class StackedBlocksFactory : public LevelContactNetworkFactory<GridType, dims>{
private:
diff --git a/src/io/hdf5-writer.hh b/src/io/hdf5-writer.hh
index 211cb8ef12b5e0a2bac29e851c4ee1e7924cf538..d0351be71cdc9e4e61c82fa32520ef81b5ba1707 100644
--- a/src/io/hdf5-writer.hh
+++ b/src/io/hdf5-writer.hh
@@ -23,7 +23,6 @@ class HDF5Writer {
HDF5Writer(HDF5::Grouplike &file,
const std::vector<Vector>& vertexCoordinates,
const std::vector<const VertexBasis* >& vertexBases,
- const std::vector<Patch>& surfaces,
const std::vector<Patch>& frictionalBoundaries,
const std::vector<ConvexPolyhedron<LocalVector>>& weakPatches)
: bodyCount_(vertexCoordinates.size()),
@@ -33,14 +32,12 @@ class HDF5Writer {
#if MY_DIM == 3
patchInfoWriters_(bodyCount_),
#endif
- surfaceWriters_(bodyCount_),
frictionalBoundaryWriters_(bodyCount_) {
for (size_t i=0; i<bodyCount_; i++) {
#if MY_DIM == 3
patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, *vertexBases[i], frictionalBoundaries[i], weakPatches[i], i);
#endif
- surfaceWriters_[i] = new SurfaceWriter<ProgramState, GridView>(file_, vertexCoordinates[i], surfaces[i], i);
frictionalBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates[i], frictionalBoundaries[i], i);
}
}
@@ -50,7 +47,6 @@ class HDF5Writer {
#if MY_DIM == 3
delete patchInfoWriters_[i];
#endif
- delete surfaceWriters_[i];
delete frictionalBoundaryWriters_[i];
}
}
@@ -65,7 +61,6 @@ class HDF5Writer {
#if MY_DIM == 3
patchInfoWriters_[i]->write(programState);
#endif
- surfaceWriters_[i]->write(programState);
frictionalBoundaryWriters_[i]->write(programState, *friction[i]);
}
}
@@ -85,7 +80,6 @@ class HDF5Writer {
#if MY_DIM == 3
std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_;
#endif
- std::vector<SurfaceWriter<ProgramState, GridView>* > surfaceWriters_;
std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionalBoundaryWriters_;
};
#endif
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index c98e933debf34bae0af1a0707a3eda16887a9035..8212c5812d2a6e871c051438b06893b8b91cf11b 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 1 # 2e-3 [m]
+smallestDiameter = 0.5 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5
diff --git a/src/multi-body-problem-data/grid/cube.cc b/src/multi-body-problem-data/grid/cube.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e3161f49ee94b46be4a0e2e0089310bf29c589b0
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cube.cc
@@ -0,0 +1,72 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "cube.hh"
+
+template <int dim>
+Cube<dim>::Cube(bool invariant = false) : invariant_(invariant) {
+ if (invariant_) {
+ nChildren_ = 1;
+ } else {
+ nChildren_ = std::pow(2, dim);
+ }
+ children_.resize(nChildren_, nullptr);
+}
+
+template <int dim>
+Cube<dim>::Cube(const Vector& A, const Vector& B, bool invariant = false) : invariant_(invariant) {
+ setCorners(A, B);
+
+ if (invariant_) {
+ nChildren_ = 1;
+ } else {
+ nChildren_ = std::pow(2, dim);
+ }
+ children_.resize(nChildren_, nullptr);
+}
+
+// generate a combination of unit basis vectors from binary number
+template <int dim>
+void Cube<dim>::parseBinary(int binary, Vector& shift) const {
+
+ shift = 0;
+ for (int d=0; d<dim; d++) {
+ // yields true, if d-th digit of binary is a 1
+ if (binary & std::pow(2, d)) {
+ // use d-th basis vector
+ shift[d] = 1;
+ }
+ }
+}
+
+// constructs child cubes and sets children_
+template <int dim>
+void Cube<dim>::split() {
+ if (invariant_) {
+ children_[0] = this;
+ } else {
+ const int nNewCubes = std::pow(2, dim);
+ newCubes.resize(nNewCubes);
+
+ Vector midPoint = A_;
+ midPoint += 1/2*(B_ - A_);
+
+ const double h = std::pow(2.0, std::round(std::log(midPoint[0]-A_[0])/std::log(2)));
+
+ newCubes[0] = new Cube(A_, midPoint, false);
+ newCubes[nNewCubes-1] = new Cube(midPoint, B_, true);
+
+ for (int i=1; i<nNewCubes-1; i++) {
+ Vector shift;
+ parseBinary(i, shift);
+ shift *= h;
+
+ newCubes[i] = new Cube(A_+shift, midPoint+shift);
+ }
+ } else {
+ children_.resize(1);
+ }
+}
+
+#include "cube_tmpl.cc"
diff --git a/src/multi-body-problem-data/grid/cube.hh b/src/multi-body-problem-data/grid/cube.hh
new file mode 100644
index 0000000000000000000000000000000000000000..56ecafe8f736e3b729a8887d2d5778061a49f838
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cube.hh
@@ -0,0 +1,49 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBE_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBE_HH
+
+
+// works in any space dimension
+template <int dim>
+class Cube {
+private:
+ using Vector = Dune::FieldVector<double, dim>;
+
+ // generate a combination of unit basis vectors from binary number
+ void parseBinary(int binary, Vector& shift) const;
+
+public:
+ Cube(bool invariant = false);
+ Cube(const Vector& A, const Vector& B, bool invariant = false);
+
+ void setCorners(const Vector& A, const Vector& B) {
+ A_ = A;
+ B_ = B;
+ }
+
+ void setParent(Cube* parent) {
+ parent_ = parent;
+ }
+
+ const std::vector<std::shared_ptr<Cube>>& children() const {
+ return children_;
+ }
+
+ // constructs child cubes and sets children_
+ void split();
+
+private:
+ using Vector = Dune::FieldVector<double, dim>;
+
+ Vector A_; // two corners that are diagonal define dim-cube in any space dimension
+ Vector B_;
+
+ const bool invariant_; // flag to set if Cube can be split()
+
+ std::shared_ptr<Cube<dim>> parent_ = nullptr;
+ std::vector<std::shared_ptr<Cube<dim>>> children_;
+ int nChildren_;
+};
+
+#endif
+
+
diff --git a/src/multi-body-problem-data/grid/cube_tmpl.cc b/src/multi-body-problem-data/grid/cube_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d7aa099e3a1f7052805726d5ff74b6bc650a4c15
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cube_tmpl.cc
@@ -0,0 +1,5 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+template class Cube<MY_DIM>;
diff --git a/src/multi-body-problem-data/grid/cubefaces.cc b/src/multi-body-problem-data/grid/cubefaces.cc
new file mode 100644
index 0000000000000000000000000000000000000000..255057c8d44aefce6b2f041be2e3ab650d858dcf
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cubefaces.cc
@@ -0,0 +1,84 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "cubefaces.hh"
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+ Vector const &c) {
+ return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+ Vector const &x) {
+ auto const minmax0 = std::minmax(v1[0], v2[0]);
+ auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+ if (minmax0.first - 1e-14 * lengthScale_ > x[0] or
+ x[0] > minmax0.second + 1e-14 * lengthScale_)
+ return false;
+ if (minmax1.first - 1e-14 * lengthScale_ > x[1] or
+ x[1] > minmax1.second + 1e-14 * lengthScale_)
+ return false;
+
+ return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+ Vector const &x) {
+ return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+CubeFaces<GridView>::CubeFaces(const GridView& gridView, const Cube& cube, double lengthScale)
+ :
+ #if MY_DIM == 3
+ lower(gridView),
+ right(gridView),
+ upper(gridView),
+ left(gridView),
+ front(gridView),
+ back(gridView),
+ #else
+ lower(gridView),
+ right(gridView),
+ upper(gridView),
+ left(gridView),
+ #endif
+ cube_(cube),
+ lengthScale_(lengthScale)
+ {
+ lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(cuboidGeometry.A, cuboidGeometry.B, in.geometry().center());
+ });
+
+ right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(cuboidGeometry.B, cuboidGeometry.C, in.geometry().center());
+ });
+
+ upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(cuboidGeometry.D, cuboidGeometry.C, in.geometry().center());
+ });
+
+ left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(cuboidGeometry.A, cuboidGeometry.D, in.geometry().center());
+ });
+
+ #if MY_DIM == 3
+ front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return isClose(cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+ });
+
+ back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return isClose(-cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+ });
+ #endif
+ }
+
+#include "cubefaces_tmpl.cc"
diff --git a/src/multi-body-problem-data/grid/cubefaces.hh b/src/multi-body-problem-data/grid/cubefaces.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4262f24f9d444585e9fde1c0f4807888eb5f0edf
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cubefaces.hh
@@ -0,0 +1,76 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBEFACES_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBEFACES_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "cube.hh"
+
+template <class GridView>
+class CubeFaces {
+private:
+ using Cube = Cube<GridView::dimensionworld>;
+
+ bool isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14 * lengthScale_;
+ }
+
+ bool isClose2(double a, double b) {
+ return std::abs(a - b) <
+ 1e-14 * lengthScale_ * lengthScale_;
+ }
+
+ template <class Vector>
+ bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
+
+ template <class Vector>
+ bool xyCollinear(Vector const &a, Vector const &b, Vector const &c);
+
+ template <class Vector>
+ bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
+
+public:
+ CubeFaces(const GridView& gridView, const Cube<GridView::dimensionworld>& cube, double lengthScale);
+
+ const BoundaryPatch<GridView>& lower() const {
+ return lower_;
+ }
+
+ const BoundaryPatch<GridView>& right() const {
+ return right_;
+ }
+
+ const BoundaryPatch<GridView>& upper() const {
+ return upper_;
+ }
+
+ const BoundaryPatch<GridView>& left() const {
+ return left_;
+ }
+
+#if MY_DIM == 3
+ const BoundaryPatch<GridView>& front() const {
+ return front_;
+ }
+ const BoundaryPatch<GridView>& back() const {
+ return back_;
+ }
+#endif
+
+private:
+ BoundaryPatch<GridView> lower_;
+ BoundaryPatch<GridView> right_;
+ BoundaryPatch<GridView> upper_;
+ BoundaryPatch<GridView> left_;
+
+#if MY_DIM == 3
+ BoundaryPatch<GridView> front_;
+ BoundaryPatch<GridView> back_;
+#endif
+
+ const Cube& cube_;
+ const double lengthScale_;
+};
+
+#endif
+
+
diff --git a/src/multi-body-problem-data/grid/cubefaces_tmpl.cc b/src/multi-body-problem-data/grid/cubefaces_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..aea2ea22b7b3b8bb8e398ca507b0a5ce6c3bd06f
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cubefaces_tmpl.cc
@@ -0,0 +1,8 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+
+template struct CubeFaces<DeformedGrid::LeafGridView>;
+template struct CubeFaces<DeformedGrid::LevelGridView>;
diff --git a/src/multi-body-problem-data/grid/cubegridconstructor.hh b/src/multi-body-problem-data/grid/cubegridconstructor.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7cb6b6a30211b4fc9196d1419f57cf73f940a54d
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cubegridconstructor.hh
@@ -0,0 +1,28 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBEGRIDCONSTRUCTOR_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBEGRIDCONSTRUCTOR_HH
+
+#include "gridconstructor.hh"
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+
+template <class GridType>
+class CubeGridConstructor : public GridConstructor {
+public:
+ using Cube = Cube<GridType::dimensionworld>;
+
+ CubeGridConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
+
+ template <class GridView>
+ void constructFaces(const GridView& gridView, CuboidGeometry const &cuboidGeometry, CubeFaces<GridView>& cubeFaces);
+
+private:
+ std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries;
+};
+
+
+#endif
+
+
diff --git a/src/multi-body-problem-data/cuboidgeometry.cc b/src/multi-body-problem-data/grid/cuboidgeometry.cc
similarity index 100%
rename from src/multi-body-problem-data/cuboidgeometry.cc
rename to src/multi-body-problem-data/grid/cuboidgeometry.cc
diff --git a/src/multi-body-problem-data/cuboidgeometry.hh b/src/multi-body-problem-data/grid/cuboidgeometry.hh
similarity index 98%
rename from src/multi-body-problem-data/cuboidgeometry.hh
rename to src/multi-body-problem-data/grid/cuboidgeometry.hh
index d189502eda9e007c7ebab67697c2e8a3eaba4eb7..834f3f4d330f92d77b498dbef42ab3194a64c02b 100644
--- a/src/multi-body-problem-data/cuboidgeometry.hh
+++ b/src/multi-body-problem-data/grid/cuboidgeometry.hh
@@ -3,7 +3,7 @@
#include <dune/common/fvector.hh>
-#include "midpoint.hh"
+#include "../midpoint.hh"
class CuboidGeometry {
public:
diff --git a/src/multi-body-problem-data/grid/gridconstructor.hh b/src/multi-body-problem-data/grid/gridconstructor.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4faa6fe1254c0745eff55850217d562b1e4ceb7b
--- /dev/null
+++ b/src/multi-body-problem-data/grid/gridconstructor.hh
@@ -0,0 +1,65 @@
+#ifndef SRC_GRIDCONSTRUCTOR_HH
+#define SRC_GRIDCONSTRUCTOR_HH
+
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "../utils/diameter.hh"
+
+template <class GridType>
+class GridConstructor {
+ private:
+ double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
+ return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
+ }
+
+ public:
+ GridConstructor() {}
+
+ std::shared_ptr<GridType> grid() {
+ return grid_;
+ }
+
+ template <class LocalVector>
+ void refine(const ConvexPolyhedron<LocalVector>& weakPatch, double smallestDiameter, double lengthScale) {
+ bool needRefine = true;
+ while (true) {
+ needRefine = false;
+ for (auto &&e : elements(grid_->leafGridView())) {
+ auto const geometry = e.geometry();
+
+ auto const weakeningRegionDistance =
+ distance(weakPatch, geometry, 1e-6 * lengthScale);
+ auto const admissibleDiameter =
+ computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter, lengthScale);
+
+ if (diameter(geometry) <= admissibleDiameter)
+ continue;
+
+ needRefine = true;
+ grid_->mark(1, e);
+ }
+ if (!needRefine)
+ break;
+
+ grid_->preAdapt();
+ grid_->adapt();
+ grid_->postAdapt();
+ }
+ }
+
+ virtual void createGrid() = 0;
+
+ private:
+ Dune::GridFactory<GridType> gridFactory_;
+ std::shared_ptr<GridType> grid_;
+};
+
+
+
+
+#endif
+
+
diff --git a/src/multi-body-problem-data/mygrids.cc b/src/multi-body-problem-data/grid/mygrids.cc
similarity index 83%
rename from src/multi-body-problem-data/mygrids.cc
rename to src/multi-body-problem-data/grid/mygrids.cc
index a754ee56ce8f229196fa25fd9e3849bc4ef1284b..3d2fafda4cdf9ec5fce0a216ec63c85be52f6ec9 100644
--- a/src/multi-body-problem-data/mygrids.cc
+++ b/src/multi-body-problem-data/grid/mygrids.cc
@@ -5,48 +5,10 @@
#include <dune/fufem/geometry/polyhedrondistance.hh>
#include "mygrids.hh"
-#include "midpoint.hh"
-#include "../utils/diameter.hh"
+#include "../midpoint.hh"
+#include "../../utils/diameter.hh"
-#if MY_DIM == 3
-SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
-#endif
-
-// back-to-front, front-to-back, front-to-back
-void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
- unsigned int W) {
-#if MY_DIM == 3
- unsigned int const U2 = U + shift_;
- unsigned int const V2 = V + shift_;
- unsigned int const W2 = W + shift_;
-
- simplices_.push_back({ U, V, W, U2 });
- simplices_.push_back({ V, V2, W2, U2 });
- simplices_.push_back({ W, W2, U2, V });
-#else
- simplices_.push_back({ U, V, W });
-#endif
-}
-
-// back-to-front, back-to-front, front-to-back
-void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
- unsigned int W) {
-#if MY_DIM == 3
- unsigned int const U2 = U + shift_;
- unsigned int const V2 = V + shift_;
- unsigned int const W2 = W + shift_;
-
- simplices_.push_back({ U, V, W, U2 });
- simplices_.push_back({ V, V2, W, U2 });
- simplices_.push_back({ V2, W, U2, W2 });
-#else
- simplices_.push_back({ U, V, W });
-#endif
-}
-
-auto SimplexManager::getSimplices() -> SimplexList const & {
- return simplices_;
-}
+#include "simplexmanager.hh"
// Fix: 3D case (still Elias code)
template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_) :
diff --git a/src/multi-body-problem-data/mygrids.hh b/src/multi-body-problem-data/grid/mygrids.hh
similarity index 81%
rename from src/multi-body-problem-data/mygrids.hh
rename to src/multi-body-problem-data/grid/mygrids.hh
index aa62f7f8dc4789e121e84208fca9ce26c768979c..f37716ebb57f3f51656f747e8fdf5696e1f4aa69 100644
--- a/src/multi-body-problem-data/mygrids.hh
+++ b/src/multi-body-problem-data/grid/mygrids.hh
@@ -27,12 +27,12 @@ template <class GridView> struct MyFaces {
bool isClose(double a, double b) {
return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale;
- };
+ }
bool isClose2(double a, double b) {
return std::abs(a - b) <
1e-14 * cuboidGeometry.lengthScale * cuboidGeometry.lengthScale;
- };
+ }
template <class Vector>
bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
@@ -44,27 +44,6 @@ template <class GridView> struct MyFaces {
bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
};
-class SimplexManager {
-public:
- using SimplexList = std::vector<std::vector<unsigned int>>;
-
-#if MY_DIM == 3
- SimplexManager(unsigned int shift);
-#endif
-
- void addFromVerticesFBB(unsigned int U, unsigned int V, unsigned int W);
- void addFromVerticesFFB(unsigned int U, unsigned int V, unsigned int W);
-
- SimplexList const &getSimplices();
-
-private:
- SimplexList simplices_;
-
-#if MY_DIM == 3
- unsigned int const shift_;
-#endif
-};
-
template <class Grid> class GridsConstructor {
public:
GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
diff --git a/src/multi-body-problem-data/mygrids_tmpl.cc b/src/multi-body-problem-data/grid/mygrids_tmpl.cc
similarity index 91%
rename from src/multi-body-problem-data/mygrids_tmpl.cc
rename to src/multi-body-problem-data/grid/mygrids_tmpl.cc
index fc8ccba75f356ab8e8cebd2041693f868ccf3fb9..32a1975979054b67089c8940a51a8762327f0927 100644
--- a/src/multi-body-problem-data/mygrids_tmpl.cc
+++ b/src/multi-body-problem-data/grid/mygrids_tmpl.cc
@@ -2,8 +2,8 @@
#error MY_DIM unset
#endif
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
#include "cuboidgeometry.hh"
template class GridsConstructor<Grid>;
diff --git a/src/multi-body-problem-data/grid/simplexmanager.cc b/src/multi-body-problem-data/grid/simplexmanager.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4e4a3aaf7e07d37b106a0c464331989f7c071dc7
--- /dev/null
+++ b/src/multi-body-problem-data/grid/simplexmanager.cc
@@ -0,0 +1,45 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "simplexmanager.hh"
+
+#if MY_DIM == 3
+SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
+#endif
+
+// back-to-front, front-to-back, front-to-back
+void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
+ unsigned int W) {
+#if MY_DIM == 3
+ unsigned int const U2 = U + shift_;
+ unsigned int const V2 = V + shift_;
+ unsigned int const W2 = W + shift_;
+
+ simplices_.push_back({ U, V, W, U2 });
+ simplices_.push_back({ V, V2, W2, U2 });
+ simplices_.push_back({ W, W2, U2, V });
+#else
+ simplices_.push_back({ U, V, W });
+#endif
+}
+
+// back-to-front, back-to-front, front-to-back
+void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
+ unsigned int W) {
+#if MY_DIM == 3
+ unsigned int const U2 = U + shift_;
+ unsigned int const V2 = V + shift_;
+ unsigned int const W2 = W + shift_;
+
+ simplices_.push_back({ U, V, W, U2 });
+ simplices_.push_back({ V, V2, W, U2 });
+ simplices_.push_back({ V2, W, U2, W2 });
+#else
+ simplices_.push_back({ U, V, W });
+#endif
+}
+
+auto SimplexManager::getSimplices() -> SimplexList const & {
+ return simplices_;
+}
diff --git a/src/multi-body-problem-data/grid/simplexmanager.hh b/src/multi-body-problem-data/grid/simplexmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..74217833e98af9c0bc7f2a571c9e5fa0391769e0
--- /dev/null
+++ b/src/multi-body-problem-data/grid/simplexmanager.hh
@@ -0,0 +1,29 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_SIMPLEXMANAGER_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_SIMPLEXMANAGER_HH
+
+#include <vector>
+
+class SimplexManager {
+public:
+ using SimplexList = std::vector<std::vector<unsigned int>>;
+
+#if MY_DIM == 3
+ SimplexManager(unsigned int shift);
+#endif
+
+ void addFromVerticesFBB(unsigned int U, unsigned int V, unsigned int W);
+ void addFromVerticesFFB(unsigned int U, unsigned int V, unsigned int W);
+
+ auto getSimplices() -> SimplexList const &;
+
+private:
+ SimplexList simplices_;
+
+#if MY_DIM == 3
+ unsigned int const shift_;
+#endif
+};
+
+#endif
+
+
diff --git a/src/multi-body-problem-data/mybody.hh b/src/multi-body-problem-data/mybody.hh
index a79e359ee7e6670fed506bd4a64c2885cad7e5e8..c0160dde33964101263ba0b70cd13ac7f4ce55d4 100644
--- a/src/multi-body-problem-data/mybody.hh
+++ b/src/multi-body-problem-data/mybody.hh
@@ -8,7 +8,7 @@
#include <dune/tectonic/bodydata.hh>
#include <dune/tectonic/gravity.hh>
-#include "cuboidgeometry.hh"
+#include "grid/cuboidgeometry.hh"
#include "segmented-function.hh"
template <int dimension> class MyBodyData : public BodyData<dimension> {
diff --git a/src/multi-body-problem-data/segmented-function.hh b/src/multi-body-problem-data/segmented-function.hh
index a75a7199559a5a47c8b8527a9a9fef13a6e842f0..866415af358cfb19939af139b70e33b4f652e92b 100644
--- a/src/multi-body-problem-data/segmented-function.hh
+++ b/src/multi-body-problem-data/segmented-function.hh
@@ -5,7 +5,7 @@
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
-#include "cuboidgeometry.hh"
+#include "grid/cuboidgeometry.hh"
class SegmentedFunction
: public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
diff --git a/src/multi-body-problem-data/weakpatch.hh b/src/multi-body-problem-data/weakpatch.hh
index 4e54105bc293a3142033c1e34b04f58ad58c9540..a1fde1cb46c1e7351e650b061ae6ca584199c9d3 100644
--- a/src/multi-body-problem-data/weakpatch.hh
+++ b/src/multi-body-problem-data/weakpatch.hh
@@ -5,7 +5,7 @@
#include <dune/fufem/geometry/convexpolyhedron.hh>
-#include "cuboidgeometry.hh"
+#include "grid/cuboidgeometry.hh"
template <class LocalVector>
void getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry, ConvexPolyhedron<LocalVector>& lowerWeakPatch, ConvexPolyhedron<LocalVector>& upperWeakPatch) {
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index f02c02aaffff5563a256630e307ed5e709159865..de02c49787df22f05ef2cc28ec4b38749533d36d 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -67,7 +67,7 @@
#include "multi-body-problem-data/bc.hh"
#include "multi-body-problem-data/mybody.hh"
-#include "multi-body-problem-data/mygrids.hh"
+#include "multi-body-problem-data/grid/mygrids.hh"
//#include "spatial-solving/solverfactory.hh"
@@ -144,8 +144,13 @@ int main(int argc, char *argv[]) {
/* Vector def(levelContactNetwork.deformedGrids()[i]->size(dims));
def = 1;
deformedGridComplex.setDeformation(def, i);*/
+ auto& levelViews = levelContactNetwork.levelViews(i);
- writeToVTK(levelContactNetwork.leafView(i), "", "deformedGrid" + std::to_string(i));
+ for (size_t j=0; j<levelViews.size(); j++) {
+ writeToVTK(*levelViews[j], "", "body_" + std::to_string(i) + "_level_" + std::to_string(j));
+ }
+
+ writeToVTK(levelContactNetwork.leafView(i), "", "body_" + std::to_string(i) + "_leaf");
}
// ----------------------------
@@ -193,7 +198,10 @@ int main(int argc, char *argv[]) {
// ------------------------
levelContactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
- /*
+ typename LevelContactNetwork::GlobalFrictionContainer globalFriction;
+ levelContactNetwork.globalFriction(globalFriction);
+ globalFriction.updateAlpha(programState.alpha);
+
using MyVertexBasis = typename Assembler::VertexBasis;
using MyCellBasis = typename Assembler::CellBasis;
std::vector<Vector> vertexCoordinates(bodyCount);
@@ -214,14 +222,16 @@ int main(int argc, char *argv[]) {
vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
+ typename LevelContactNetwork::BoundaryPatches frictionBoundaries;
+ levelContactNetwork.boundaryPatches("friction", frictionBoundaries);
auto dataWriter =
writeData ? std::make_unique<
HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
*dataFile, vertexCoordinates, vertexBases,
- surfaces, frictionalBoundaries, weakPatches)
+ frictionBoundaries, weakPatches)
: nullptr;
-
+/*
const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "body");
IterationRegister iterationCount;
@@ -261,34 +271,58 @@ int main(int argc, char *argv[]) {
}
};
report(true);
-
-
+*/
// -------------------
// Set up TNNMG solver
// -------------------
auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
- const auto& totalDirichletNodes = levelContactNetwork.totalDirichletNodes();
+
+ Dune::BitSetVector<dims> totalDirichletNodes;
+ levelContactNetwork.totalNodes("dirichlet", totalDirichletNodes);
using NonlinearFactory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, totalDirichletNodes);
- using MyUpdater = Updaters<RateUpdater<Vector, Matrix, typename LevelContactNetwork::Body::Function, dims>,
+ using BoundaryFunctions = typename LevelContactNetwork::BoundaryFunctions;
+ using BoundaryNodes = typename LevelContactNetwork::BoundaryNodes;
+ 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"));
- MyUpdater current(
- initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
- velocityDirichletFunctions, dirichletNodes, matrices,
- programState.u, programState.v, programState.a),
- initStateUpdater<ScalarVector, Vector>(
+
+ BoundaryFunctions velocityDirichletFunctions;
+ levelContactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+ BoundaryNodes dirichletNodes;
+ levelContactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+ typename LevelContactNetwork::BoundaryPatchNodes frictionNodes;
+ levelContactNetwork.boundaryPatchNodes("friction", frictionNodes);
+
+ Updaters current(
+ initRateUpdater(
+ parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunctions,
+ dirichletNodes,
+ levelContactNetwork.matrices(),
+ programState.u,
+ programState.v,
+ programState.a),
+ initStateUpdater<ScalarVector, Vector, typename LevelContactNetwork::BoundaryPatchNodes>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
- programState.alpha, frictionalNodes,
- L, V0));
+ programState.alpha,
+ frictionNodes,
+ L,
+ V0));
auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
- auto const mustRefine = [&](MyUpdater &coarseUpdater,
- MyUpdater &fineUpdater) {
+
+ typename LevelContactNetwork::StateEnergyNorms stateEnergyNorms;
+ levelContactNetwork.stateEnergyNorms(stateEnergyNorms);
+
+ auto const mustRefine = [&](Updaters &coarseUpdater,
+ Updaters &fineUpdater) {
std::vector<ScalarVector> coarseAlpha;
coarseUpdater.state_->extractAlpha(coarseAlpha);
@@ -307,11 +341,14 @@ int main(int argc, char *argv[]) {
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
- AdaptiveTimeStepper<NonlinearFactory, MyUpdater, EnergyNorm<ScalarMatrix, ScalarVector>>
+ typename LevelContactNetwork::ExternalForces externalForces;
+ levelContactNetwork.externalForces(externalForces);
+
+ AdaptiveTimeStepper<NonlinearFactory, Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
adaptiveTimeStepper(factory, parset, globalFriction, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
-
+/*
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index 03c7084ae23e1d8d4a14c76f78504531efe42000..b850b5cc829734e50342d1e812254778f74ad97b 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -11,7 +11,7 @@ vtk.write = false
[problem]
finalTime = 1000 # [s]
-bodyCount = 3
+bodyCount = 2
[body]
bulkModulus = 0.5e5 # [Pa]
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index c3e90842cc708f2375e613b71e53f70eed241a8d..6b776feb982debaca743bd3bb6598d4a8ebdd8ca 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -36,7 +36,7 @@ void FixedPointIterationCounter::operator+=(
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
Factory &factory, Dune::ParameterTree const &parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, const std::vector<const ErrorNorm* >& errorNorms)
+ GlobalFrictionContainer& globalFriction, const std::vector<const ErrorNorm* >& errorNorms)
: factory_(factory),
step_(factory_.getStep()),
parset_(parset),
@@ -81,9 +81,8 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
relativeVelocities(velocityIterates, v_rel);
// contribution from nonlinearity
- for (size_t i=0; i<alpha.size(); i++) {
- globalFriction_[i]->updateAlpha(alpha[i]);
- }
+ globalFriction_.updateAlpha(alpha);
+
std::vector<Matrix> matrices(velocityMatrices.size());
std::vector<Vector> rhs(velocityRHSs.size());
@@ -91,8 +90,11 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
matrices[i] = velocityMatrices[i];
rhs[i] = velocityRHSs[i];
- globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
- globalFriction_[i]->addGradient(v_rel[i], rhs[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];
}
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index bae4624f1dc8f8c04e7de03ac4165ac2bf438aae..688a80b516c3a8dfe4e34764a9a48688477118a7 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -10,6 +10,8 @@
#include <dune/contact/assemblers/nbodyassembler.hh>
+#include "../data-structures/globalfrictioncontainer.hh"
+
struct FixedPointIterationCounter {
size_t iterations = 0;
size_t multigridIterations = 0;
@@ -28,15 +30,17 @@ class FixedPointIterator {
using Nonlinearity = typename Factory::Nonlinearity;
// using Nonlinearity = typename ConvexProblem::NonlinearityType;
+ using DeformedGrid = typename Factory::DeformedGrid;
- using DeformedGrid = typename Factory::DeformedGrid;
+public:
+ using GlobalFrictionContainer = GlobalFrictionContainer<Nonlinearity, 2>;
private:
void relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const;
public:
FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
+ GlobalFrictionContainer& globalFriction,
const std::vector<const ErrorNorm* >& errorNorms);
FixedPointIterationCounter run(Updaters updaters,
@@ -48,7 +52,7 @@ class FixedPointIterator {
Factory& factory_;
std::shared_ptr<typename Factory::Step> step_;
Dune::ParameterTree const &parset_;
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
+ GlobalFrictionContainer& globalFriction_;
size_t fixedPointMaxIterations_;
double fixedPointTolerance_;
diff --git a/src/spatial-solving/fixedpointiterator_tmpl.cc b/src/spatial-solving/fixedpointiterator_tmpl.cc
index 97dfe8341850194d93d6cba4dca335331c3342d6..eec3461083148f6a5b8b372dcd7abe3cd655e61a 100644
--- a/src/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/src/spatial-solving/fixedpointiterator_tmpl.cc
@@ -15,16 +15,21 @@
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/myblockproblem.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"
-using Function = Dune::VirtualFunction<double, double>;
+
+
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+
using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 2d83a0c13aa36ce8c2c5ea006cec0cdc7a4b7730..0bba4b88001d23557b8608aaef3f7f40f02fd42f 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -26,6 +26,11 @@ SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactory
Solver::QUIET),
transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
+ /* baseSolverStep_.obstacles_ = using HasObstacle = Dune::BitSetVector<BlockSize>;
+ using Obstacle = BoxConstraint<Field, BlockSize>;
+ const HasObstacle* hasObstacle_;
+ const std::vector<Obstacle>* obstacles_;*/
+
multigridStep_ = std::make_shared<Step>();
// tnnmg iteration step
multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
@@ -34,7 +39,7 @@ SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactory
multigridStep_->setSmoother(&presmoother_, &postsmoother_);
multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
-
+ multigridStep_->setVerbosity(NumProc::QUIET);
// create the transfer operators
const int nCouplings = nBodyAssembler_.nCouplings();
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 9e982c90b3367e1cb5a1e84410d2e4367145a0b0..5e87ea27d7920dbab7b044a485c5a7d4f8f159e7 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -6,8 +6,10 @@
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
+
//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
@@ -50,7 +52,8 @@ class SolverFactory {
private:
const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler_;
- ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
+ //ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
+ BlockGSStep<MatrixType, VectorType> baseSolverStep_;
EnergyNorm<MatrixType, VectorType> baseEnergyNorm_;
LoopSolver<VectorType> baseSolver_;
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4109f1b97ceead8006cffc009c8ee1cbd2817826
--- /dev/null
+++ b/src/tests/CMakeLists.txt
@@ -0,0 +1 @@
+dune_add_test(SOURCES globalfrictioncontainertest.cc)
diff --git a/src/tests/globalfrictioncontainertest.cc b/src/tests/globalfrictioncontainertest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7439c584284f180e22313407c4d142fb0ea9c877
--- /dev/null
+++ b/src/tests/globalfrictioncontainertest.cc
@@ -0,0 +1,120 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <vector>
+
+#include <exception>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+
+#include "../data-structures/globalfrictioncontainer.hh"
+
+struct Body {
+ size_t id;
+ size_t localAlpha;
+
+ void updateAlpha(size_t lAlpha) {
+ localAlpha = lAlpha;
+ }
+};
+
+int main(int argc, char *argv[]) {
+ try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::cout << "--------------------------" << std::endl;
+ std::cout << "-- GlobalFriction Test: --" << std::endl;
+ std::cout << "--------------------------" << std::endl << std::endl;
+
+ // have to be <10
+ size_t n = 5;
+ size_t m = 3;
+ size_t l = 2;
+
+ std::vector<std::vector<int>> ref;
+ ref.resize(n);
+ for (size_t i=0; i<n; i++) {
+ ref[i].resize(m);
+ for (size_t j=0; j<m; j++) {
+ ref[i][j] = 10*i + j;
+ }
+ }
+
+ // resize() test
+ bool resizeTest = true;
+ {
+ GlobalFrictionContainer<Body,1> globalFriction1;
+ globalFriction1.resize({n});
+ resizeTest = resizeTest & (globalFriction1.size() == n);
+
+ GlobalFrictionContainer<Body,2> globalFriction2;
+ globalFriction2.resize({n});
+ resizeTest = resizeTest & (globalFriction2.size() == n);
+
+ for (size_t i=0; i<globalFriction2.size(); i++) {
+ resizeTest = resizeTest & (globalFriction2[i].size() == 0);
+ }
+ }
+
+ // operator[] test
+ bool operatorTest = true;
+ {
+ GlobalFrictionContainer<Body,2> globalFriction;
+ globalFriction.resize({n, m});
+ operatorTest = operatorTest & (globalFriction.size() == n);
+
+ for (size_t i=0; i<n; i++) {
+ for (size_t j=0; j<m; j++) {
+ operatorTest = operatorTest & (globalFriction[i][j] == nullptr);
+ }
+ }
+ }
+
+ // updateAlpha() test
+ bool updateTest = true;
+ {
+ GlobalFrictionContainer<Body,3> globalFriction;
+ globalFriction.resize({n,m,l});
+
+ for (size_t i=0; i<n; i++) {
+ for (size_t j=0; j<m; j++) {
+ for (size_t k=0; k<l; k++) {
+ globalFriction[i][j][k] = std::make_shared<Body>();
+ globalFriction[i][j][k]->id = 10*i + j;
+ }
+ }
+ }
+
+ globalFriction.updateAlpha(ref);
+
+ for (size_t i=0; i<n; i++) {
+ for (size_t j=0; j<m; j++) {
+ for (size_t k=0; k<l; k++) {
+ updateTest = updateTest & (globalFriction[i][j][k]->id == globalFriction[i][j][k]->localAlpha);
+ }
+ }
+ }
+
+ }
+
+ std::cout << "resize(): " << resizeTest << std::endl;
+ std::cout << "operator[]: " << operatorTest << std::endl;
+ std::cout << "updateAlpha(): " << updateTest << std::endl;
+ std::cout << "-------------------------" << std::endl << std::endl;
+
+ bool passed = resizeTest & operatorTest & updateTest;
+
+ std::cout << "Overall the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+ return passed ? 0 : 1;
+
+ } catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+ } catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+ }
+}
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index 1c36768ee8fb049be3a357b57d053c5a58e6a61c..cde9e81d81cb6d445e4acdc8aec166be90c56d45 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -20,9 +20,9 @@ void IterationRegister::reset() {
template <class Factory, class Updaters, class ErrorNorm>
AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters ¤t,
+ GlobalFrictionContainer& globalFriction, Updaters ¤t,
double relativeTime, double relativeTau,
- std::vector<std::function<void(double, Vector &)>>& externalForces,
+ ExternalForces& externalForces,
const std::vector<const ErrorNorm* >& errorNorms,
std::function<bool(Updaters &, Updaters &)> mustRefine)
: relativeTime_(relativeTime),
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index f7897a2906dbf22d53cc5bf334abb329cc437ddc..7b251fade0bb1709a209bacae861c696593a6050 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -24,16 +24,19 @@ class AdaptiveTimeStepper {
using Vector = typename Factory::Vector;
//using ConvexProblem = typename Factory::ConvexProblem;
- using Nonlinearity = typename Factory::Nonlinearity;
+ //using Nonlinearity = typename Factory::Nonlinearity;
using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
+ using GlobalFrictionContainer = typename MyCoupledTimeStepper::GlobalFrictionContainer;
+ using ExternalForces = typename MyCoupledTimeStepper::ExternalForces;
+
public:
AdaptiveTimeStepper(Factory &factory, Dune::ParameterTree const &parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
+ GlobalFrictionContainer& globalFriction,
Updaters ¤t, double relativeTime,
double relativeTau,
- std::vector<std::function<void(double, Vector &)>>& externalForces,
+ ExternalForces& externalForces,
const std::vector<const ErrorNorm* >& errorNorms,
std::function<bool(Updaters &, Updaters &)> mustRefine);
@@ -50,10 +53,10 @@ class AdaptiveTimeStepper {
double finalTime_;
Factory &factory_;
Dune::ParameterTree const &parset_;
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
+ GlobalFrictionContainer& globalFriction_;
Updaters ¤t_;
UpdatersWithCount R1_;
- std::vector<std::function<void(double, Vector &)>>& externalForces_;
+ ExternalForces& externalForces_;
std::function<bool(Updaters &, Updaters &)> mustRefine_;
const std::vector<const ErrorNorm* >& errorNorms_;
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index 97c383edff112abcd9c782cc6d50bc2860cf3e33..4f11f38799fb9caa92e811342a2ac15da900a5e5 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -13,15 +13,18 @@
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/myblockproblem.hh>
+#include "../data-structures/levelcontactnetwork.hh"
#include "../spatial-solving/solverfactory.hh"
#include "rate/rateupdater.hh"
#include "state/stateupdater.hh"
#include "updaters.hh"
-using Function = Dune::VirtualFunction<double, double>;
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+
using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 74e66fc80594a18a69db83932038372c538882b0..e06ce73233e419f737b8455c354b3fe3ad5e0c61 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -7,9 +7,9 @@
template <class Factory, class Updaters, class ErrorNorm>
CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
double finalTime, Factory &factory, Dune::ParameterTree const &parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters updaters,
+ GlobalFrictionContainer& globalFriction, Updaters updaters,
const std::vector<const ErrorNorm* >& errorNorms,
- std::vector<std::function<void(double, Vector &)>>& externalForces)
+ ExternalForces& externalForces)
: finalTime_(finalTime),
factory_(factory),
parset_(parset),
@@ -28,7 +28,7 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
auto const newRelativeTime = relativeTime + relativeTau;
std::vector<Vector> ell(externalForces_.size());
for (size_t i=0; i<externalForces_.size(); i++) {
- externalForces_[i](newRelativeTime, ell[i]);
+ (*externalForces_[i])(newRelativeTime, ell[i]);
}
std::vector<Matrix> velocityMatrix;
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 97c0128243e64a9431c7a80ea9ce9dc1eab1342c..81a53c6fd5bd286dee6951a4eb6a64fbae9fb923 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -12,14 +12,17 @@ template <class Factory, class Updaters, class ErrorNorm>
class CoupledTimeStepper {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
- using Nonlinearity = typename Factory::Nonlinearity;
+ //using Nonlinearity = typename Factory::Nonlinearity;
+public:
+ using GlobalFrictionContainer = typename FixedPointIterator<Factory, Updaters, ErrorNorm>::GlobalFrictionContainer;
+ using ExternalForces = std::vector<const std::function<void(double, Vector &)>*>;
public:
CoupledTimeStepper(double finalTime, Factory &factory,
Dune::ParameterTree const &parset,
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
+ GlobalFrictionContainer& globalFriction,
Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
- std::vector<std::function<void(double, Vector &)>>& externalForces);
+ ExternalForces& externalForces);
FixedPointIterationCounter step(double relativeTime, double relativeTau);
@@ -27,9 +30,9 @@ class CoupledTimeStepper {
double finalTime_;
Factory &factory_;
Dune::ParameterTree const &parset_;
- std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
+ GlobalFrictionContainer& globalFriction_;
Updaters updaters_;
- std::vector<std::function<void(double, Vector &)>>& externalForces_;
+ ExternalForces& externalForces_;
const std::vector<const ErrorNorm* >& errorNorms_;
};
#endif
diff --git a/src/time-stepping/coupledtimestepper_tmpl.cc b/src/time-stepping/coupledtimestepper_tmpl.cc
index a2bd2dfa43a144219f95c7793de0520f85fea1c5..e16e9df9c948a02153eb818a75334410cd1fd88b 100644
--- a/src/time-stepping/coupledtimestepper_tmpl.cc
+++ b/src/time-stepping/coupledtimestepper_tmpl.cc
@@ -13,16 +13,19 @@
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/myblockproblem.hh>
+#include "../data-structures/levelcontactnetwork.hh"
#include "../spatial-solving/solverfactory.hh"
#include "rate/rateupdater.hh"
#include "state/stateupdater.hh"
#include "updaters.hh"
-using Function = Dune::VirtualFunction<double, double>;
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+
using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
diff --git a/src/time-stepping/rate.cc b/src/time-stepping/rate.cc
index b94a2a353cef977101784289a2a88633119a65eb..5ff412aac0f72b4b9b81315e2b8e72716ebe8435 100644
--- a/src/time-stepping/rate.cc
+++ b/src/time-stepping/rate.cc
@@ -6,23 +6,24 @@
#include "rate/backward_euler.hh"
#include "rate/newmark.hh"
-template <class Vector, class Matrix, class Function, int dimension>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
-initRateUpdater(Config::scheme scheme,
- const std::vector<const Function*>& velocityDirichletFunctions,
- const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto initRateUpdater(Config::scheme scheme,
+ const BoundaryFunctions& velocityDirichletFunctions,
+ const BoundaryNodes& velocityDirichletNodes,
const Matrices<Matrix,2>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
- const std::vector<Vector>& a_initial) {
+ const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
switch (scheme) {
case Config::Newmark:
- return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
+ return std::make_shared<Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(
matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
velocityDirichletFunctions);
case Config::BackwardEuler:
return std::make_shared<
- BackwardEuler<Vector, Matrix, Function, dimension>>(
+ BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(
matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
velocityDirichletFunctions);
default:
diff --git a/src/time-stepping/rate.hh b/src/time-stepping/rate.hh
index 8a3cc285779d0918bc93c25066eb4a924039e99e..4497a4e1e4c598182402d8d00a46c00c7910cbcf 100644
--- a/src/time-stepping/rate.hh
+++ b/src/time-stepping/rate.hh
@@ -6,13 +6,13 @@
#include "../data-structures/enums.hh"
#include "rate/rateupdater.hh"
-template <class Vector, class Matrix, class Function, int dimension>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
-initRateUpdater(Config::scheme scheme,
- const std::vector<const Function* >& velocityDirichletFunctions,
- const std::vector<std::vector<Dune::BitSetVector<dimension>>>& velocityDirichletNodes,
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto initRateUpdater(Config::scheme scheme,
+ const BoundaryFunctions& velocityDirichletFunctions,
+ const BoundaryNodes& velocityDirichletNodes,
const Matrices<Matrix,2>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
- const std::vector<Vector>& a_initial);
+ const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>;
#endif
diff --git a/src/time-stepping/rate/backward_euler.cc b/src/time-stepping/rate/backward_euler.cc
index c4845782846180f98dabf1283415bef265ec6ad3..b7a915f9c1e38b9252543e72af4de1384198b94f 100644
--- a/src/time-stepping/rate/backward_euler.cc
+++ b/src/time-stepping/rate/backward_euler.cc
@@ -3,102 +3,121 @@
#include "backward_euler.hh"
-template <class Vector, class Matrix, class Function, size_t dim>
-BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
- const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
- const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
- const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<const Function*>& _dirichletFunctions)
- : RateUpdater<Vector, Matrix, Function, dim>(
- _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
- _dirichletFunctions) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
- double _tau,
- double relativeTime,
- std::vector<Vector>& rhs, std::vector<Vector>& iterate,
- std::vector<Matrix>& AM) {
- for (size_t i=0; i<this->u.size(); i++) {
- this->dirichletFunctions[i]->evaluate(relativeTime, this->dirichletValue);
- this->tau = _tau;
-
- /* We start out with the formulation
-
- M a + C v + A u = ell
-
- Backward Euler means
-
- a1 = 1.0/tau ( v1 - v0 )
- u1 = tau v1 + u0
-
- in summary, we get at time t=1
-
- M [1.0/tau ( v1 - v0 )] + C v1
- + A [tau v1 + u0] = ell
-
- or
-
- 1.0/tau M v1 + C v1 + tau A v1
- = [1.0/tau M + C + tau A] v1
- = ell + 1.0/tau M v0 - A u0
- */
-
- // set up LHS (for fixed tau, we'd only really have to do this once)
- Matrix& LHS = AM[i];
- {
- Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
- this->matrices.elasticity[i]->M());
- indices.import(*this->matrices.elasticity[i]);
- indices.import(*this->matrices.mass[i]);
- indices.import(*this->matrices.damping[i]);
- indices.exportIdx(LHS);
- }
- LHS = 0.0;
- Dune::MatrixVector::addProduct(LHS, 1.0 / this->tau, *this->matrices.mass[i]);
- Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
- Dune::MatrixVector::addProduct(LHS, this->tau, *this->matrices.elasticity[i]);
-
- // set up RHS
- {
- Vector& rhss = rhs[i];
- rhss = ell[i];
- Dune::MatrixVector::addProduct(rhss, 1.0 / this->tau, *this->matrices.mass[i],
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::BackwardEuler(
+ const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const BoundaryNodes& _dirichletNodes,
+ const BoundaryFunctions& _dirichletFunctions) :
+ RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+ _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
+ _dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
+ const std::vector<Vector>& ell,
+ double _tau,
+ double relativeTime,
+ std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+ std::vector<Matrix>& AM) {
+
+ for (size_t i=0; i<this->u.size(); i++) {
+ this->tau = _tau;
+
+ /* We start out with the formulation
+
+ M a + C v + A u = ell
+
+ Backward Euler means
+
+ a1 = 1.0/tau ( v1 - v0 )
+ u1 = tau v1 + u0
+
+ in summary, we get at time t=1
+
+ M [1.0/tau ( v1 - v0 )] + C v1
+ + A [tau v1 + u0] = ell
+
+ or
+
+ 1.0/tau M v1 + C v1 + tau A v1
+ = [1.0/tau M + C + tau A] v1
+ = ell + 1.0/tau M v0 - A u0
+ */
+
+ // set up LHS (for fixed tau, we'd only really have to do this once)
+ Matrix& LHS = AM[i];
+ {
+ Dune::MatrixIndexSet indices(
+ this->matrices.elasticity[i]->N(),
+ this->matrices.elasticity[i]->M());
+ indices.import(*this->matrices.elasticity[i]);
+ indices.import(*this->matrices.mass[i]);
+ indices.import(*this->matrices.damping[i]);
+ indices.exportIdx(LHS);
+ }
+ LHS = 0.0;
+ Dune::MatrixVector::addProduct(LHS, 1.0 / this->tau, *this->matrices.mass[i]);
+ Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+ Dune::MatrixVector::addProduct(LHS, this->tau, *this->matrices.elasticity[i]);
+
+ // set up RHS
+ {
+ Vector& rhss = rhs[i];
+ rhss = ell[i];
+ Dune::MatrixVector::addProduct(rhss, 1.0 / this->tau, *this->matrices.mass[i],
this->v_o[i]);
- Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+ Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+ }
+
}
iterate = this->v_o;
- const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
- for (size_t k = 0; k < dirichletNodess.size(); ++k)
- for (size_t j = 0; j < dim; ++j)
- if (dirichletNodess[k][j])
- iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
- }
+ for (size_t i=0; i<iterate.size(); i++) {
+ auto& bodyIterate = iterate[i];
+
+ const auto& bodyDirichletNodes = this->dirichletNodes[i];
+ const auto& bodyDirichletFunctions = this->dirichletFunctions[i];
+ for (size_t bc=0; bc<bodyDirichletNodes.size(); ++bc) {
+ const auto& bcDirichletNodes = *bodyDirichletNodes[bc];
+ size_t dim = bcDirichletNodes[0].size();
+
+ double dirichletValue;
+ bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
+
+ for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
+ for (size_t j=0; j<dim; ++j) {
+ if (bcDirichletNodes[k][j]) {
+ iterate[k][j] = dirichletValue;
+ }
+ }
+ }
+ }
+ }
}
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
- const std::vector<Vector>& iterate) {
- this->postProcessCalled = true;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::postProcess(const std::vector<Vector>& iterate) {
+ this->postProcessCalled = true;
- this->v = iterate;
+ this->v = iterate;
- this->u = this->u_o;
+ this->u = this->u_o;
- for (size_t i=0; i<this->u.size(); i++) {
- Dune::MatrixVector::addProduct(this->u[i], this->tau, this->v[i]);
+ for (size_t i=0; i<this->u.size(); i++) {
+ Dune::MatrixVector::addProduct(this->u[i], this->tau, this->v[i]);
- Vector& ai = this->a[i];
- ai = this->v[i];
- ai -= this->v_o[i];
- ai /= this->tau;
- }
+ Vector& ai = this->a[i];
+ ai = this->v[i];
+ ai -= this->v_o[i];
+ ai /= this->tau;
+ }
}
-template <class Vector, class Matrix, class Function, size_t dim>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>>
-BackwardEuler<Vector, Matrix, Function, dim>::clone() const {
- return std::make_shared<BackwardEuler<Vector, Matrix, Function, dim>>(*this);
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::clone() const
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
+ return std::make_shared<BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(*this);
}
diff --git a/src/time-stepping/rate/backward_euler.hh b/src/time-stepping/rate/backward_euler.hh
index 26e5831d26a184a51d17af729837225294dca122..bb48a2550bb64672b4ab0092b91764e46346477b 100644
--- a/src/time-stepping/rate/backward_euler.hh
+++ b/src/time-stepping/rate/backward_euler.hh
@@ -1,22 +1,27 @@
#ifndef SRC_TIME_STEPPING_RATE_BACKWARD_EULER_HH
#define SRC_TIME_STEPPING_RATE_BACKWARD_EULER_HH
-template <class Vector, class Matrix, class Function, size_t dim>
-class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class BackwardEuler : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
public:
- BackwardEuler(const Matrices<Matrix,2> &_matrices, const std::vector<Vector> &_u_initial,
- const std::vector<Vector> &_v_initial, const std::vector<Vector> &_a_initial,
- const std::vector<Dune::BitSetVector<dim> > &_dirichletNodes,
- const std::vector<const Function*> &_dirichletFunctions);
+ BackwardEuler(
+ const Matrices<Matrix,2> &_matrices,
+ const std::vector<Vector> &_u_initial,
+ const std::vector<Vector> &_v_initial,
+ const std::vector<Vector> &_a_initial,
+ const BoundaryNodes& _dirichletNodes,
+ const BoundaryFunctions& _dirichletFunctions);
- void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
- std::vector<Matrix>&) override;
+ virtual void setup(
+ const std::vector<Vector>&,
+ double,
+ double,
+ std::vector<Vector>&,
+ std::vector<Vector>&,
+ std::vector<Matrix>&) override;
- void postProcess(const std::vector<Vector>&) override;
+ virtual void postProcess(const std::vector<Vector>&) override;
- std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
- const override;
-
-private:
+ virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
};
#endif
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index b3fb46699fb58f959a6e0ef891ffc0cdca3bae06..2961e7e8885a92681c35fe0b8f48aa7d76d2cee4 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -3,108 +3,127 @@
#include "newmark.hh"
-template <class Vector, class Matrix, class Function, size_t dim>
-Newmark<Vector, Matrix, Function, dim>::Newmark(
- const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
- const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
- const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<const Function*>& _dirichletFunctions)
- : RateUpdater<Vector, Matrix, Function, dim>(
- _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
- _dirichletFunctions) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
- double _tau,
- double relativeTime,
- std::vector<Vector>& rhs, std::vector<Vector>& iterate,
- std::vector<Matrix>& AM) {
- for (size_t i=0; i<this->u.size(); i++) {
- this->dirichletFunctions[i]->evaluate(relativeTime, this->dirichletValue);
- this->tau = _tau;
-
- /* We start out with the formulation
-
- M a + C v + A u = ell
-
- Newmark means
-
- a1 = 2/tau ( v1 - v0 ) - a0
- u1 = tau/2 ( v1 + v0 ) + u0
-
- in summary, we get at time t=1
-
- M [2/tau ( u1 - u0 ) - a0] + C v1
- + A [tau/2 ( v1 + v0 ) + u0] = ell
-
- or
-
- 2/tau M v1 + C v1 + tau/2 A v1
- = [2/tau M + C + tau/2 A] v1
- = ell + 2/tau M v0 + M a0
- - tau/2 A v0 - A u0
- */
-
- // set up LHS (for fixed tau, we'd only really have to do this once)
- Matrix& LHS = AM[i];
- {
- Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
- this->matrices.elasticity[i]->M());
- indices.import(*this->matrices.elasticity[i]);
- indices.import(*this->matrices.mass[i]);
- indices.import(*this->matrices.damping[i]);
- indices.exportIdx(LHS);
- }
- LHS = 0.0;
- Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
- Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
- Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
-
- // set up RHS
- {
- Vector& rhss = rhs[i];
- rhss = ell[i];
- Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::Newmark(
+ const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const BoundaryNodes& _dirichletNodes,
+ const BoundaryFunctions& _dirichletFunctions) :
+ RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+ _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
+ _dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
+ const std::vector<Vector>& ell,
+ double _tau,
+ double relativeTime,
+ std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+ std::vector<Matrix>& AM) {
+
+ for (size_t i=0; i<this->u.size(); i++) {
+ this->tau = _tau;
+
+ /* We start out with the formulation
+
+ M a + C v + A u = ell
+
+ Newmark means
+
+ a1 = 2/tau ( v1 - v0 ) - a0
+ u1 = tau/2 ( v1 + v0 ) + u0
+
+ in summary, we get at time t=1
+
+ M [2/tau ( u1 - u0 ) - a0] + C v1
+ + A [tau/2 ( v1 + v0 ) + u0] = ell
+
+ or
+
+ 2/tau M v1 + C v1 + tau/2 A v1
+ = [2/tau M + C + tau/2 A] v1
+ = ell + 2/tau M v0 + M a0
+ - tau/2 A v0 - A u0
+ */
+
+ // set up LHS (for fixed tau, we'd only really have to do this once)
+ Matrix& LHS = AM[i];
+ {
+ Dune::MatrixIndexSet indices(
+ this->matrices.elasticity[i]->N(),
+ this->matrices.elasticity[i]->M());
+ indices.import(*this->matrices.elasticity[i]);
+ indices.import(*this->matrices.mass[i]);
+ indices.import(*this->matrices.damping[i]);
+ indices.exportIdx(LHS);
+ }
+ LHS = 0.0;
+ Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
+ Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+ Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
+
+ // set up RHS
+ {
+ Vector& rhss = rhs[i];
+ rhss = ell[i];
+ Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
+ this->v_o[i]);
+ Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
+ Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
this->v_o[i]);
- Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
- Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
- this->v_o[i]);
- Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+ Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+ }
}
iterate = this->v_o;
- const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
- for (size_t k = 0; k < dirichletNodess.size(); ++k)
- for (size_t j = 0; j < dim; ++j)
- if (dirichletNodess[k][j])
- iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
- }
+ for (size_t i=0; i<iterate.size(); i++) {
+ auto& bodyIterate = iterate[i];
+
+ const auto& bodyDirichletNodes = this->dirichletNodes[i];
+ const auto& bodyDirichletFunctions = this->dirichletFunctions[i];
+ for (size_t bc=0; bc<bodyDirichletNodes.size(); ++bc) {
+ const auto& bcDirichletNodes = *bodyDirichletNodes[bc];
+ size_t dim = bcDirichletNodes[0].size();
+
+ double dirichletValue;
+ bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
+
+ for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
+ for (size_t j=0; j<dim; ++j) {
+ if (bcDirichletNodes[k][j]) {
+ iterate[k][j] = dirichletValue;
+ }
+ }
+ }
+ }
+ }
}
-template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::postProcess(const std::vector<Vector>& iterate) {
- this->postProcessCalled = true;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::postProcess(const std::vector<Vector>& iterate) {
+ this->postProcessCalled = true;
- this->v = iterate;
+ this->v = iterate;
- // u1 = tau/2 ( v1 + v0 ) + u0
- this->u = this->u_o;
+ // u1 = tau/2 ( v1 + v0 ) + u0
+ this->u = this->u_o;
- for (size_t i=0; i<this->u.size(); i++) {
- Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
- Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
+ for (size_t i=0; i<this->u.size(); i++) {
+ Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
+ Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
- // a1 = 2/tau ( v1 - v0 ) - a0
- this->a[i] = 0.0;
- Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
- Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
- Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
- }
+ // a1 = 2/tau ( v1 - v0 ) - a0
+ this->a[i] = 0.0;
+ Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
+ Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
+ Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
+ }
}
-template <class Vector, class Matrix, class Function, size_t dim>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>>
-Newmark<Vector, Matrix, Function, dim>::clone() const {
- return std::make_shared<Newmark<Vector, Matrix, Function, dim>>(*this);
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::clone() const
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
+ return std::make_shared<Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(*this);
}
diff --git a/src/time-stepping/rate/newmark.hh b/src/time-stepping/rate/newmark.hh
index 9041e01dc066dde3a03b9285930aaf595598af95..4aef6caf6b2c24ac36754c647350e1f68f641c60 100644
--- a/src/time-stepping/rate/newmark.hh
+++ b/src/time-stepping/rate/newmark.hh
@@ -1,20 +1,27 @@
#ifndef SRC_TIME_STEPPING_RATE_NEWMARK_HH
#define SRC_TIME_STEPPING_RATE_NEWMARK_HH
-template <class Vector, class Matrix, class Function, size_t dim>
-class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class Newmark : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
public:
- Newmark(const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
- const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
- const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<const Function*>& _dirichletFunctions);
+ Newmark(
+ const Matrices<Matrix,2>& _matrices,
+ const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial,
+ const std::vector<Vector>& _a_initial,
+ const BoundaryNodes& _dirichletNodes,
+ const BoundaryFunctions& _dirichletFunctions);
- void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
- std::vector<Matrix>&) override;
+ virtual void setup(
+ const std::vector<Vector>&,
+ double,
+ double,
+ std::vector<Vector>&,
+ std::vector<Vector>&,
+ std::vector<Matrix>&) override;
- void postProcess(const std::vector<Vector>&) override;
+ virtual void postProcess(const std::vector<Vector>&) override;
- std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
- const override;
+ virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
};
#endif
diff --git a/src/time-stepping/rate/rateupdater.cc b/src/time-stepping/rate/rateupdater.cc
index fb2e22a2d78caa630774df09d9d1aa79cd1b52c2..aa29295b7b8e6d8ac8c209e5d484beb624ad3531 100644
--- a/src/time-stepping/rate/rateupdater.cc
+++ b/src/time-stepping/rate/rateupdater.cc
@@ -4,12 +4,12 @@
#include "rateupdater.hh"
-template <class Vector, class Matrix, class Function, size_t dim>
-RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::RateUpdater(
const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
- const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<const Function*>& _dirichletFunctions)
+ const BoundaryNodes& _dirichletNodes,
+ const BoundaryFunctions& _dirichletFunctions)
: matrices(_matrices),
u(_u_initial),
v(_v_initial),
@@ -17,41 +17,41 @@ RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
dirichletNodes(_dirichletNodes),
dirichletFunctions(_dirichletFunctions) {}
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
- u_o = u;
- v_o = v;
- a_o = a;
- postProcessCalled = false;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::nextTimeStep() {
+ u_o = u;
+ v_o = v;
+ a_o = a;
+ postProcessCalled = false;
}
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(std::vector<Vector>& displacements) const {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractDisplacement(std::vector<Vector>& displacements) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
- displacements = u;
+ displacements = u;
}
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(std::vector<Vector>& velocity) const {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractVelocity(std::vector<Vector>& velocity) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
- velocity = v;
+ velocity = v;
}
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
- oldVelocity = v_o;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
+ oldVelocity = v_o;
}
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(std::vector<Vector>& acceleration) const {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractAcceleration(std::vector<Vector>& acceleration) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
- acceleration = a;
+ acceleration = a;
}
#include "backward_euler.cc"
diff --git a/src/time-stepping/rate/rateupdater.hh b/src/time-stepping/rate/rateupdater.hh
index a9819ff5894f14673cc944451b774676426ce227..5c06af42f058a98912699ecbf3b99b4b772499f0 100644
--- a/src/time-stepping/rate/rateupdater.hh
+++ b/src/time-stepping/rate/rateupdater.hh
@@ -7,18 +7,26 @@
#include "../../data-structures/matrices.hh"
-template <class Vector, class Matrix, class Function, size_t dim>
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
class RateUpdater {
protected:
- RateUpdater(const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
- const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
- const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
- const std::vector<const Function*>& _dirichletFunctions);
+ RateUpdater(
+ const Matrices<Matrix,2>& _matrices,
+ const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial,
+ const std::vector<Vector>& _a_initial,
+ const BoundaryNodes& _dirichletNodes,
+ const BoundaryFunctions& _dirichletFunctions);
public:
- void nextTimeStep();
- void virtual setup(const std::vector<Vector>& ell, double _tau, double relativeTime,
- std::vector<Vector>& rhs, std::vector<Vector>& iterate, std::vector<Matrix>& AB) = 0;
+ void nextTimeStep();
+ void virtual setup(
+ const std::vector<Vector>& ell,
+ double _tau,
+ double relativeTime,
+ std::vector<Vector>& rhs,
+ std::vector<Vector>& iterate,
+ std::vector<Matrix>& AB) = 0;
void virtual postProcess(const std::vector<Vector>& iterate) = 0;
void extractDisplacement(std::vector<Vector>& displacements) const;
@@ -26,15 +34,13 @@ class RateUpdater {
void extractOldVelocity(std::vector<Vector>& velocity) const;
void extractAcceleration(std::vector<Vector>& acceleration) const;
- std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone()
- const = 0;
+ std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> virtual clone() const = 0;
protected:
const Matrices<Matrix,2>& matrices;
std::vector<Vector> u, v, a;
- const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
- const std::vector<const Function*>& dirichletFunctions;
- double dirichletValue;
+ const BoundaryNodes& dirichletNodes;
+ const BoundaryFunctions& dirichletFunctions;
std::vector<Vector> u_o, v_o, a_o;
double tau;
diff --git a/src/time-stepping/rate/rateupdater_tmpl.cc b/src/time-stepping/rate/rateupdater_tmpl.cc
index 0af5138618c73ba1fba451fe77a6b8cd6f555472..273c261198849795ad1fd0d9f39b316b80f256b2 100644
--- a/src/time-stepping/rate/rateupdater_tmpl.cc
+++ b/src/time-stepping/rate/rateupdater_tmpl.cc
@@ -2,12 +2,14 @@
#error MY_DIM unset
#endif
-#include <dune/common/function.hh>
-
+#include "../../explicitgrid.hh"
#include "../../explicitvectors.hh"
-using Function = Dune::VirtualFunction<double, double>;
+#include "../../data-structures/levelcontactnetwork.hh"
+
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
-template class RateUpdater<Vector, Matrix, Function, MY_DIM>;
-template class Newmark<Vector, Matrix, Function, MY_DIM>;
-template class BackwardEuler<Vector, Matrix, Function, MY_DIM>;
+template class RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index b71b59eb3e85b3e0255f0bff47bcbd5a92adba37..0d1271ed0e11cf146900a568c435c6716ce521da 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -1,15 +1,16 @@
+#include "../data-structures/levelcontactnetwork_tmpl.cc"
#include "../explicitvectors.hh"
-#include <dune/common/function.hh>
+using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryFunctions;
+using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryNodes;
-using Function = Dune::VirtualFunction<double, double>;
-
-template std::shared_ptr<RateUpdater<Vector, Matrix, Function, MY_DIM>>
-initRateUpdater<Vector, Matrix, Function, MY_DIM>(
+template
+auto initRateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
Config::scheme scheme,
- const std::vector<const Function* >& velocityDirichletFunctions,
- const std::vector<Dune::BitSetVector<MY_DIM>>& velocityDirichletNodes,
+ const BoundaryFunctions& velocityDirichletFunctions,
+ const BoundaryNodes& velocityDirichletNodes,
const Matrices<Matrix, 2>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
- const std::vector<Vector>& a_initial);
+ const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>;
diff --git a/src/time-stepping/state.cc b/src/time-stepping/state.cc
index b8e4349eae5dead583fede83ed929845b018e186..02233a4468bb254cf8d41117c30c2e02dab31403 100644
--- a/src/time-stepping/state.cc
+++ b/src/time-stepping/state.cc
@@ -6,16 +6,22 @@
#include "state/ageinglawstateupdater.cc"
#include "state/sliplawstateupdater.cc"
-template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0) {
+template <class ScalarVector, class Vector, class BoundaryPatchNodes>
+auto initStateUpdater(
+ Config::stateModel model,
+ const std::vector<ScalarVector>& alpha_initial,
+ const BoundaryPatchNodes& frictionNodes,
+ const std::vector<double>& L,
+ const std::vector<double>& V0)
+-> std::shared_ptr<StateUpdater<ScalarVector, Vector>> {
+
switch (model) {
case Config::AgeingLaw:
return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionalNodes, L, V0);
+ alpha_initial, frictionNodes, L, V0);
case Config::SlipLaw:
return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionalNodes, L, V0);
+ alpha_initial, frictionNodes, L, V0);
default:
assert(false);
}
diff --git a/src/time-stepping/state.hh b/src/time-stepping/state.hh
index 3c70b2db5d0b04263d194e6dc9b29a580962cfa7..ca8352ad6bb73c35aa2425a846f5e373e40edcd3 100644
--- a/src/time-stepping/state.hh
+++ b/src/time-stepping/state.hh
@@ -10,9 +10,12 @@
#include "state/sliplawstateupdater.hh"
#include "state/stateupdater.hh"
-template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
- Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
-
+template <class ScalarVector, class Vector, class BoundaryPatchNodes>
+auto initStateUpdater(
+ Config::stateModel model,
+ const std::vector<ScalarVector>& alpha_initial,
+ const BoundaryPatchNodes& frictionNodes,
+ const std::vector<double>& L,
+ const std::vector<double>& V0)
+-> std::shared_ptr<StateUpdater<ScalarVector, Vector>>;
#endif
diff --git a/src/time-stepping/state/ageinglawstateupdater.cc b/src/time-stepping/state/ageinglawstateupdater.cc
index 1d7bb0ee570dd34c657684fcaac8f99a71a7b79a..94a7075ceb5f1fcb7d9965168a48e7a45581de51 100644
--- a/src/time-stepping/state/ageinglawstateupdater.cc
+++ b/src/time-stepping/state/ageinglawstateupdater.cc
@@ -2,14 +2,39 @@
#include "ageinglawstateupdater.hh"
#include "../../utils/tobool.hh"
+#include "../../utils/debugutils.hh"
template <class ScalarVector, class Vector>
+template <class BoundaryPatchNodes>
AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
- const std::vector<ScalarVector>& _alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& _nodes,
- const std::vector<double>& _L,
- const std::vector<double>& _V0)
- : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
+ 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;
+ }
+ }
+ }
+ }
+}
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
@@ -25,7 +50,10 @@ void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
Compute [ 1-\exp(c*x) ] / x under the assumption that x is
non-negative
*/
-double liftSingularity(double c, double x) {
+auto liftSingularity(
+ double c,
+ double x) {
+
if (x <= 0)
return -c;
else
@@ -33,34 +61,36 @@ double liftSingularity(double c, double x) {
}
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) +
+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));
+ }
}
- }
}
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
- std::vector<ScalarVector>& _alpha) {
- _alpha = alpha;
+ std::vector<ScalarVector>& _alpha) {
+
+ _alpha = alpha;
}
template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-AgeingLawStateUpdater<ScalarVector, Vector>::clone() const {
+auto AgeingLawStateUpdater<ScalarVector, Vector>::clone() const
+ -> std::shared_ptr<StateUpdater<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 4b8ec5588167b412a76d65a5477b988d10e3a12b..0c9ab450341bb2288dd41bdddd6d996144e24c90 100644
--- a/src/time-stepping/state/ageinglawstateupdater.hh
+++ b/src/time-stepping/state/ageinglawstateupdater.hh
@@ -6,21 +6,24 @@
template <class ScalarVector, class Vector>
class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
public:
- AgeingLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
- const std::vector<double>& _V0);
+ template <class BoundaryPatchNodes>
+ AgeingLawStateUpdater(
+ const std::vector<ScalarVector>& _alpha_initial,
+ const BoundaryPatchNodes& _nodes,
+ const std::vector<double>& _L,
+ const std::vector<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;
- std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+ auto clone() const -> std::shared_ptr<StateUpdater<ScalarVector, Vector>> override;
private:
std::vector<ScalarVector> alpha_o;
std::vector<ScalarVector> alpha;
- const std::vector<Dune::BitSetVector<1>>& nodes;
+ std::vector<Dune::BitSetVector<1>> nodes;
const std::vector<double>& L;
const std::vector<double>& V0;
double tau;
diff --git a/src/time-stepping/state/sliplawstateupdater.cc b/src/time-stepping/state/sliplawstateupdater.cc
index 196b02b09be73a622e9983441d09932be1ff8b7c..51ae920bd2ed5d1d51e51dc4f7b5bd9ee39153cc 100644
--- a/src/time-stepping/state/sliplawstateupdater.cc
+++ b/src/time-stepping/state/sliplawstateupdater.cc
@@ -4,52 +4,78 @@
#include "../../utils/tobool.hh"
template <class ScalarVector, class Vector>
+template <class BoundaryPatchNodes>
SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
- const std::vector<ScalarVector>& _alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
- const std::vector<double>& _V0)
- : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
+ 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];
+ auto& _bodyNodes = _nodes[i];
+ if (_bodyNodes.size() > 0) {
+ bodyNodes.resize(_bodyNodes[0]->size(), false);
+ }
+
+ 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;
+ }
+ }
+ }
+ }
+}
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;
+ tau = _tau;
}
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::solve(
- const std::vector<Vector>& velocity_field) {
+ 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 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;
+ 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]) +
+ 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);
- }
- }
+ }
+ }
}
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
std::vector<ScalarVector>& _alpha) {
+
_alpha = alpha;
}
template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-SlipLawStateUpdater<ScalarVector, Vector>::clone() const {
+auto SlipLawStateUpdater<ScalarVector, Vector>::clone() const
+-> std::shared_ptr<StateUpdater<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 760044e83f12fefe7fce144c3c4378ff647177da..0b3f7ef39170b6db0f572bd7429b94ab9090dc87 100644
--- a/src/time-stepping/state/sliplawstateupdater.hh
+++ b/src/time-stepping/state/sliplawstateupdater.hh
@@ -6,21 +6,24 @@
template <class ScalarVector, class Vector>
class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
public:
- SlipLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
- const std::vector<double>& _V0);
+ template <class BoundaryPatchNodes>
+ SlipLawStateUpdater(
+ const std::vector<ScalarVector>& _alpha_initial,
+ const BoundaryPatchNodes& _nodes,
+ const std::vector<double>& _L,
+ const std::vector<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;
- std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+ auto clone() const -> std::shared_ptr<StateUpdater<ScalarVector, Vector>> override;
private:
std::vector<ScalarVector> alpha_o;
std::vector<ScalarVector> alpha;
- const std::vector<Dune::BitSetVector<1>>& nodes;
+ std::vector<Dune::BitSetVector<1>> nodes;
const std::vector<double>& L;
const std::vector<double>& V0;
double tau;
diff --git a/src/time-stepping/state_tmpl.cc b/src/time-stepping/state_tmpl.cc
index 718150c7af0332286c7fe489863caa35330a7c18..5046cf258356874450b8fb4ad5de5d58ed162795 100644
--- a/src/time-stepping/state_tmpl.cc
+++ b/src/time-stepping/state_tmpl.cc
@@ -1,6 +1,19 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../data-structures/levelcontactnetwork.hh"
+
+using BoundaryPatchNodes = typename LevelContactNetwork<Grid, MY_DIM>::LevelBoundaryPatchNodes;
-template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-initStateUpdater<ScalarVector, Vector>(
- Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
- const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
+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>>;
diff --git a/src/utils/almostequal.hh b/src/utils/almostequal.hh
new file mode 100644
index 0000000000000000000000000000000000000000..448c95fef53598bcad3eb0612299f2a7d31b8678
--- /dev/null
+++ b/src/utils/almostequal.hh
@@ -0,0 +1,11 @@
+#ifndef SRC_ALMOSTEQUAL_HH
+#define SRC_ALMOSTEQUAL_HH
+
+#include <math.h>
+
+template <typename ctype>
+typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type almost_equal(ctype x, ctype y, int ulp) {
+ return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp || std::abs(x-y) < std::numeric_limits<ctype>::min();
+}
+
+#endif
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index 75b5c8e8f60484be1d2a3abd826ae646e07e7228..4671ff68de50d8c51cf674d89f968e9097548f56 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -73,7 +73,6 @@ using namespace std;
std::cout << std::endl;
}
- /*
template <class T=Dune::FieldMatrix<double,1,1>>
void print(const Dune::Matrix<T>& mat, std::string message){
std::cout << message << std::endl;
@@ -95,7 +94,7 @@ using namespace std;
std::cout << std::endl << std::endl;
}
- void print(const std::vector<Dune::FieldVector<double,1>>& x, std::string message){
+ /* void print(const std::vector<Dune::FieldVector<double,1>>& x, std::string message){
std::cout << message << std::endl;
for (size_t i=0; i<x.size(); i++) {
std::cout << x[i] << " ";
@@ -156,7 +155,7 @@ using namespace std;
std::cout.rdbuf( lBufferOld );
}
-/*
+
template <class BasisType, typename ctype=double>
void printBasisDofLocation(const BasisType& basis) {
typedef typename BasisType::GridView GridView;
@@ -218,5 +217,5 @@ using namespace std;
}
std::cout << std::endl;
}
- } */
+ }
#endif