From 8e48508f58de63a8209ef4bc74585a247357d8a7 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 21 Jun 2019 18:03:28 +0200
Subject: [PATCH] .
---
src/CMakeLists.txt | 7 +-
src/data-structures/body.cc | 300 ++++++++------
src/data-structures/body.hh | 162 +++++---
src/data-structures/body_tmpl.cc | 5 +-
src/data-structures/contactnetwork.cc | 391 +++++++++++++++++-
src/data-structures/contactnetwork.hh | 153 ++++++-
src/data-structures/contactnetwork_tmpl.cc | 7 +-
src/data-structures/levelcontactnetwork.cc | 342 ++++-----------
src/data-structures/levelcontactnetwork.hh | 214 ++++------
.../levelcontactnetwork_tmpl.cc | 8 +-
src/data-structures/program_state.hh | 22 +-
src/explicitvectors.hh | 10 +-
src/factories/contactnetworkfactory.hh | 36 +-
src/factories/levelcontactnetworkfactory.hh | 8 +-
src/factories/stackedblocksfactory.cc | 1 -
src/factories/threeblocksfactory.cc | 204 +++++++++
src/factories/threeblocksfactory.hh | 84 ++++
src/factories/threeblocksfactory_tmpl.cc | 8 +
.../grid/cuboidgeometry.cc | 131 +++---
.../grid/cuboidgeometry.hh | 80 ++--
.../grid/cuboidgeometry_tmpl.cc | 8 +
src/multi-body-problem-data/grid/mygrids.cc | 30 +-
src/multi-body-problem-data/grid/mygrids.hh | 10 +-
.../grid/mygrids_tmpl.cc | 4 +-
src/multi-body-problem-data/weakpatch.hh | 67 ---
src/multi-body-problem.cc | 177 +++++---
.../contact/dualmortarcoupling.hh | 4 +-
src/spatial-solving/fixedpointiterator.cc | 14 +-
src/spatial-solving/fixedpointiterator.hh | 8 +-
.../fixedpointiterator_tmpl.cc | 12 +-
.../preconditioners/supportpatchfactory.hh | 173 ++++----
src/tests/CMakeLists.txt | 1 +
src/tests/supportpatchfactorytest.cc | 144 +++++++
src/time-stepping/adaptivetimestepper.cc | 20 +-
src/time-stepping/adaptivetimestepper.hh | 8 +-
src/time-stepping/adaptivetimestepper_tmpl.cc | 14 +-
src/time-stepping/coupledtimestepper.cc | 12 +-
src/time-stepping/coupledtimestepper.hh | 16 +-
src/time-stepping/coupledtimestepper_tmpl.cc | 12 +-
src/time-stepping/rate/rateupdater_tmpl.cc | 6 +-
src/time-stepping/rate_tmpl.cc | 6 +-
src/utils/debugutils.hh | 17 +-
42 files changed, 1953 insertions(+), 983 deletions(-)
create mode 100644 src/factories/threeblocksfactory.cc
create mode 100644 src/factories/threeblocksfactory.hh
create mode 100644 src/factories/threeblocksfactory_tmpl.cc
create mode 100644 src/multi-body-problem-data/grid/cuboidgeometry_tmpl.cc
delete mode 100644 src/multi-body-problem-data/weakpatch.hh
create mode 100644 src/tests/supportpatchfactorytest.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0cfcc708..6bf43fe7 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -28,11 +28,12 @@ set(MSW_SOURCE_FILES
assemblers.cc
#nodalweights.cc
data-structures/body.cc
- #data-structures/contactnetwork.cc
- data-structures/enumparser.cc
data-structures/levelcontactnetwork.cc
+ data-structures/contactnetwork.cc
+ data-structures/enumparser.cc
#factories/cantorfactory.cc
- factories/stackedblocksfactory.cc
+ factories/threeblocksfactory.cc
+ #factories/stackedblocksfactory.cc
io/vtk.cc
io/hdf5/frictionalboundary-writer.cc
io/hdf5/iteration-writer.cc
diff --git a/src/data-structures/body.cc b/src/data-structures/body.cc
index 39220bd5..75d92f7d 100644
--- a/src/data-structures/body.cc
+++ b/src/data-structures/body.cc
@@ -2,33 +2,31 @@
#include "config.h"
#endif
-
+#include <dune/common/hybridutilities.hh>
#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) :
+// -------------------------------
+// --- LeafBody Implementation ---
+// -------------------------------
+template <class GridTEMPLATE, class VectorType>
+LeafBody<GridTEMPLATE, VectorType>::LeafBody(
+ const std::shared_ptr<BodyData<dim>>& bodyData,
+ const std::shared_ptr<GridTEMPLATE>& hostGrid) :
bodyData_(bodyData),
- grid_(grid),
- deformation_(std::make_shared<DeformationFunction>(grid_->leafGridView())),
- deformedGrid_(std::make_shared<DeformedGridType>(*grid_, *deformation_)),
- leafView_(deformedGrid_->leafGridView()),
- leafVertexCount_(leafView_.size(dims)),
- assembler_(std::make_shared<Assembler>(leafView_)),
+ hostGrid_(hostGrid),
+ deformation_(std::make_shared<DeformationFunction>(hostGrid_->leafGridView())),
+ grid_(std::make_shared<GridType>(*hostGrid_, *deformation_)),
+ gridView_(grid_->leafGridView()),
+ assembler_(std::make_shared<Assembler>(gridView_)),
matrices_()
{
- leafBoundaryConditions_.resize(0);
- levelBoundaryConditions_.resize(grid_->maxLevel()+1);
- for (size_t i=0; i<levelBoundaryConditions_.size(); i++) {
- levelBoundaryConditions_[i].resize(0);
- }
+ boundaryConditions_.resize(0);
externalForce_ = [&](double relativeTime, VectorType &ell) {
// Problem formulation: right-hand side
- std::vector<std::shared_ptr<LeafBoundaryCondition>> leafNeumannConditions;
- leafBoundaryConditions("neumann", leafNeumannConditions);
+ std::vector<std::shared_ptr<BoundaryCondition>> leafNeumannConditions;
+ boundaryConditions("neumann", leafNeumannConditions);
ell.resize(gravityFunctional_.size());
ell = gravityFunctional_;
@@ -44,8 +42,8 @@ Body<GridTEMPLATE, VectorTEMPLATE, dims>::Body(
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assemble() {
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::assemble() {
// assemble matrices
assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
@@ -58,204 +56,274 @@ void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assemble() {
// setter and getter
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::data() const
--> const std::shared_ptr<BodyData<dims>>& {
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::data() const
+-> const std::shared_ptr<BodyData<dim>>& {
return bodyData_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::setDeformation(const VectorType& def) {
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::setDeformation(const VectorType& def) {
deformation_->setDeformation(def);
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::deformation() const
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::deformation() const
-> DeformationFunction& {
return *deformation_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::deformedGrid() const
--> const DeformedGridType* {
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::grid() const
+-> std::shared_ptr<GridType> {
- return deformedGrid_.get();
+ return grid_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafView() const
--> const DeformedLeafGridView& {
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::gridView() const
+-> const GridView& {
- return leafView_;
+ return gridView_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafVertexCount() const
--> int {
-
- return leafVertexCount_;
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::nVertices() const
+-> size_t {
+ return gridView_.size(dim);
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembler() const
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::assembler() const
-> const std::shared_ptr<Assembler>& {
return assembler_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::matrices() const
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::matrices() const
-> const Matrices<Matrix,1>& {
return matrices_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::externalForce() const
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::externalForce() const
-> const ExternalForce& {
return externalForce_;
}
-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 VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
+ boundaryConditions_.emplace_back(boundaryCondition);
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafBoundaryConditions(
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryConditions(
const std::string& tag,
- std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const {
+ std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {
selectedConditions.resize(0);
- for (size_t i=0; i<leafBoundaryConditions_.size(); i++) {
- if (leafBoundaryConditions_[i]->tag() == tag)
- selectedConditions.emplace_back(leafBoundaryConditions_[i]);
+ for (size_t i=0; i<boundaryConditions_.size(); i++) {
+ if (boundaryConditions_[i]->tag() == tag)
+ selectedConditions.emplace_back(boundaryConditions_[i]);
}
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafBoundaryConditions() const
--> const std::vector<std::shared_ptr<LeafBoundaryCondition>>& {
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::boundaryConditions() const
+-> const std::vector<std::shared_ptr<BoundaryCondition>>& {
- return leafBoundaryConditions_;
+ return boundaryConditions_;
}
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryPatchNodes(
+ const std::string& tag,
+ BoundaryPatchNodes& nodes) const {
+
+ std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(tag, _boundaryConditions);
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(
- std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition,
- int level) {
+ nodes.resize(_boundaryConditions.size());
- levelBoundaryConditions_[level].push_back(levelBoundaryCondition);
+ 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>::levelBoundaryConditions(
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryNodes(
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]);
- }
+ BoundaryNodes& nodes) const {
+
+ std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(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>::levelBoundaryConditions(
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryFunctions(
const std::string& tag,
- int level,
- std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const {
+ BoundaryFunctions& functions) 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]);
+ std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(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 VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryPatches(
+ const std::string& tag,
+ BoundaryPatches& patches) const {
+
+ std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(tag, _boundaryConditions);
+
+ patches.resize(_boundaryConditions.size());
+
+ for (size_t i=0; i<patches.size(); i++) {
+ patches[i] = _boundaryConditions[i]->boundaryPatch().get();
}
}
-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];
+// ---------------------------
+// --- Body Implementation ---
+// ---------------------------
+
+// setter and getter
+
+template <class GridView>
+auto Body<GridView>::data() const
+-> const std::shared_ptr<BodyData<dim>>& {
+
+ return bodyData_;
+}
+
+
+template <class GridView>
+auto Body<GridView>::grid() const
+-> std::shared_ptr<GridType> {
+
+ return grid_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions() const
--> const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& {
+template <class GridView>
+auto Body<GridView>::gridView() const
+-> const GridView& {
+
+ return gridView_;
+}
+
+template <class GridView>
+auto Body<GridView>::nVertices() const
+-> size_t {
+ return gridView_.size(dim);
+}
+
+
+template <class GridView>
+void Body<GridView>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
+
+ boundaryConditions_.push_back(boundaryCondition);
+}
+
+
+template <class GridView>
+void Body<GridView>::boundaryConditions(
+ const std::string& tag,
+ std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {
+
+ selectedConditions.resize(0);
+ for (size_t i=0; i<boundaryConditions_.size(); i++) {
+ if (boundaryConditions_[i]->tag() == tag)
+ selectedConditions.push_back(boundaryConditions_[i]);
+ }
+}
+
+
+template <class GridView>
+auto Body<GridView>::boundaryConditions() const
+-> const std::vector<std::shared_ptr<BoundaryCondition>>& {
- return levelBoundaryConditions_;
+ return boundaryConditions_;
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryPatchNodes(
+template <class GridView>
+void Body<GridView>::boundaryPatchNodes(
const std::string& tag,
BoundaryPatchNodes& nodes) const {
- std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
- this->leafBoundaryConditions(tag, boundaryConditions);
+ std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(tag, _boundaryConditions);
- nodes.resize(boundaryConditions.size());
+ nodes.resize(_boundaryConditions.size());
for (size_t i=0; i<nodes.size(); i++) {
- nodes[i] = boundaryConditions[i]->boundaryPatch()->getVertices();
+ nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
}
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryNodes(
+template <class GridView>
+void Body<GridView>::boundaryNodes(
const std::string& tag,
BoundaryNodes& nodes) const {
- std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
- this->leafBoundaryConditions(tag, boundaryConditions);
+ std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(tag, _boundaryConditions);
- nodes.resize(boundaryConditions.size());
+ nodes.resize(_boundaryConditions.size());
for (size_t i=0; i<nodes.size(); i++) {
- nodes[i] = boundaryConditions[i]->boundaryNodes().get();
+ nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
}
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryFunctions(
+template <class GridView>
+void Body<GridView>::boundaryFunctions(
const std::string& tag,
BoundaryFunctions& functions) const {
- std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
- this->leafBoundaryConditions(tag, boundaryConditions);
+ std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(tag, _boundaryConditions);
- functions.resize(boundaryConditions.size());
+ functions.resize(_boundaryConditions.size());
for (size_t i=0; i<functions.size(); i++) {
- functions[i] = boundaryConditions[i]->boundaryFunction().get();
+ functions[i] = _boundaryConditions[i]->boundaryFunction().get();
}
}
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryPatches(
+template <class GridView>
+void Body<GridView>::boundaryPatches(
const std::string& tag,
BoundaryPatches& patches) const {
- std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
- this->leafBoundaryConditions(tag, boundaryConditions);
+ std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+ this->boundaryConditions(tag, _boundaryConditions);
- patches.resize(boundaryConditions.size());
+ patches.resize(_boundaryConditions.size());
for (size_t i=0; i<patches.size(); i++) {
- patches[i] = boundaryConditions[i]->boundaryPatch().get();
+ patches[i] = _boundaryConditions[i]->boundaryPatch().get();
}
}
diff --git a/src/data-structures/body.hh b/src/data-structures/body.hh
index ed71b586..6b5463cd 100644
--- a/src/data-structures/body.hh
+++ b/src/data-structures/body.hh
@@ -2,6 +2,7 @@
#define SRC_DATA_STRUCTURES_BODY_HH
#include <functional>
+#include <type_traits>
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
@@ -30,51 +31,41 @@
#include "enums.hh"
#include "matrices.hh"
-
-
-template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
-class Body {
+template <class HostGridType, class VectorType>
+class LeafBody {
public:
- using GridType = GridTEMPLATE;
- using VectorType = VectorTEMPLATE;
+ static const int dim = HostGridType::dimensionworld;
- 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 HostGridType::LeafGridView, VectorType>;
+ using GridType = Dune::GeometryGrid<HostGridType, DeformationFunction>;
+ using GridView = typename GridType::LeafGridView;
- using Assembler = MyAssembler<DeformedLeafGridView, dims>;
+ using Function = Dune::VirtualFunction<double, double>;
+ using ExternalForce = std::function<void(double, VectorType &)>;
+
+ using Assembler = MyAssembler<GridView, dim>;
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 Function = Dune::VirtualFunction<double, double>;
- using ExternalForce = std::function<void(double, VectorType &)>;
-
- using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
- using LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
+ using BoundaryCondition = BoundaryCondition<GridView, dim>;
using BoundaryFunctions = std::vector<const Function* >;
- using BoundaryNodes = std::vector<const Dune::BitSetVector<dims>* >;
+ using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
- using BoundaryPatches = std::vector<const typename LeafBoundaryCondition::BoundaryPatch* >;
+ using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;
using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
private:
- std::shared_ptr<Body<GridType, VectorType, dims>> parent_;
- std::vector<std::shared_ptr<Body<GridType, VectorType, dims>>> children_;
+ std::shared_ptr<BodyData<dim>> bodyData_;
- std::shared_ptr<BodyData<dims>> bodyData_;
-
- std::shared_ptr<GridType> grid_;
+ std::shared_ptr<HostGridType> hostGrid_;
std::shared_ptr<DeformationFunction> deformation_;
- std::shared_ptr<DeformedGridType> deformedGrid_;
-
- DeformedLeafGridView leafView_;
- int leafVertexCount_;
+ std::shared_ptr<GridType> grid_;
+ GridView gridView_;
std::shared_ptr<Assembler> assembler_;
Matrices<Matrix,1> matrices_;
@@ -84,47 +75,122 @@ class Body {
VectorType gravityFunctional_;
// 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
+ std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;
public:
- Body(
- const std::shared_ptr<BodyData<dims>>& bodyData,
- const std::shared_ptr<GridType>& grid);
+ LeafBody(
+ const std::shared_ptr<BodyData<dim>>& bodyData,
+ const std::shared_ptr<HostGridType>& hostGrid);
void assemble();
// setter and getter
- auto data() const -> const std::shared_ptr<BodyData<dims>>&;
+ auto data() const -> const std::shared_ptr<BodyData<dim>>&;
void setDeformation(const VectorType& def);
auto deformation() const -> DeformationFunction&;
- auto deformedGrid() const -> const DeformedGridType*;
+ auto grid() const -> std::shared_ptr<GridType>;
- auto leafView() const -> const DeformedLeafGridView&;
- auto leafVertexCount() const -> int;
+ auto gridView() const -> const GridView&;
+ auto nVertices() const -> size_t;
auto assembler() const -> const std::shared_ptr<Assembler>&;
auto matrices() const -> const Matrices<Matrix,1>&;
auto externalForce() const -> const ExternalForce&;
- void addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition);
- void leafBoundaryConditions(
+ void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
+ void boundaryConditions(
const std::string& tag,
- std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const;
- auto leafBoundaryConditions() const -> const std::vector<std::shared_ptr<LeafBoundaryCondition>>&;
+ std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
+ auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;
- void addBoundaryCondition(std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition, int level);
- void levelBoundaryConditions(
+ void boundaryPatchNodes(
+ const std::string& tag,
+ BoundaryPatchNodes& nodes) const;
+ void boundaryNodes(
+ const std::string& tag,
+ BoundaryNodes& nodes) const;
+ void boundaryFunctions(
const std::string& tag,
- std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const;
- void levelBoundaryConditions(
+ BoundaryFunctions& functions) const;
+ void boundaryPatches(
+ const std::string& tag,
+ BoundaryPatches& patches) const;
+};
+
+template <class GridView>
+class Body {
+public:
+ enum {dim = GridView::dimensionworld};
+
+ using GridType = typename GridView::Grid;
+
+ using Function = Dune::VirtualFunction<double, double>;
+
+ using BoundaryCondition = BoundaryCondition<GridView, dim>;
+
+ using BoundaryFunctions = std::vector<const Function* >;
+ using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
+ using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+ using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;
+
+private:
+ std::shared_ptr<BodyData<dim>> bodyData_;
+
+ std::shared_ptr<GridType> grid_;
+ const size_t level_;
+ GridView gridView_;
+
+ // boundary conditions
+ std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;
+
+public:
+ template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LeafGridView>::value, int> = 0>
+ Body(
+ const std::shared_ptr<BodyData<dim>>& bodyData,
+ const std::shared_ptr<GridType>& grid) :
+ bodyData_(bodyData),
+ grid_(grid),
+ level_(grid_->maxLevel()),
+ gridView_(grid_->leafGridView()) {
+
+ boundaryConditions_.resize(0);
+ }
+
+ template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
+ Body(
+ const std::shared_ptr<BodyData<dim>>& bodyData,
+ const std::shared_ptr<GridType>& grid,
+ const size_t level) :
+ bodyData_(bodyData),
+ grid_(grid),
+ level_(level),
+ gridView_(grid_->levelGridView(level_)) {
+
+ boundaryConditions_.resize(0);
+ }
+
+
+ // setter and getter
+ auto data() const -> const std::shared_ptr<BodyData<dim>>&;
+ auto grid() const -> std::shared_ptr<GridType>;
+
+ //template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
+ auto level() const
+ -> size_t {
+
+ return level_;
+ }
+
+ auto gridView() const -> const GridView&;
+ auto nVertices() const -> size_t;
+
+ void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
+ void boundaryConditions(
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>>>&;
+ std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
+ auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;
void boundaryPatchNodes(
const std::string& tag,
diff --git a/src/data-structures/body_tmpl.cc b/src/data-structures/body_tmpl.cc
index b35e1cfa..d2a5dab1 100644
--- a/src/data-structures/body_tmpl.cc
+++ b/src/data-structures/body_tmpl.cc
@@ -5,4 +5,7 @@
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
-template class Body<Grid, Vector, MY_DIM>;
+template class LeafBody<Grid, Vector>;
+
+template class Body<DefLeafGridView>;
+template class Body<DefLevelGridView>;
diff --git a/src/data-structures/contactnetwork.cc b/src/data-structures/contactnetwork.cc
index ac118a1e..81d9f8a9 100644
--- a/src/data-structures/contactnetwork.cc
+++ b/src/data-structures/contactnetwork.cc
@@ -2,30 +2,397 @@
#include "config.h"
#endif
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/tectonic/globalratestatefriction.hh>
+#include <dune/tectonic/frictionpotential.hh>
#include "contactnetwork.hh"
-template <class GridType, int dim>
-void ContactNetwork<GridType, dim>::assemble() {
+#include "../assemblers.hh"
+#include "../utils/tobool.hh"
+
+
+template <class HostGridType, class VectorType>
+ContactNetwork<HostGridType, VectorType>::ContactNetwork(
+ size_t nBodies,
+ size_t nCouplings,
+ size_t nLevels) :
+ levelContactNetworks_(nLevels),
+ bodies_(nBodies),
+ couplings_(nCouplings),
+ frictionBoundaries_(nCouplings),
+ stateEnergyNorms_(nBodies),
+ nBodyAssembler_(nBodies, nCouplings)
+{}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::addLevel(std::shared_ptr<LevelContactNetwork> levelContactNetwork) {
+ size_t level = levelContactNetwork->level();
+ if (level >= levelContactNetworks_.size()) {
+ levelContactNetworks_.resize(level);
+ }
+
+ levelContactNetworks_[level] = levelContactNetwork;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::addLevel(const std::vector<size_t>& bodyLevels, const size_t level) {
+ assert(bodyLevels.size() == nBodies());
+
+ if (level >= levelContactNetworks_.size()) {
+ levelContactNetworks_.resize(level+1);
+ }
+
+ auto& levelContactNetwork = levelContactNetworks_[level];
+ levelContactNetwork = std::make_shared<LevelContactNetwork>(nBodies(), nCouplings(), level);
+
+ std::vector<std::shared_ptr<Body>> bodies(nBodies());
+
+ for (size_t id=0; id<nBodies(); id++) {
+ const auto& body = this->body(id);
+ bodies[id] = std::make_shared<Body>(body->data(), body->grid(), std::min((size_t) body->grid()->maxLevel(), bodyLevels[id]));
+ }
+ levelContactNetwork->setBodies(bodies);
+ levelContactNetwork->setCouplings(this->couplings_);
+ levelContactNetwork->build();
+}
+
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::assemble() {
+ std::vector<const GridType*> grids(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ grids[i] = bodies_[i]->grid().get();
+ bodies_[i]->assemble();
+ frictionBoundaries_[i] = std::make_unique<BoundaryPatch>(bodies_[i]->gridView(), false);
+ }
+
+ // set up dune-contact nBodyAssembler
+ nBodyAssembler_.setGrids(grids);
+
+ for (size_t i=0; i<nCouplings(); i++) {
+ nBodyAssembler_.setCoupling(*couplings_[i], i);
+ }
+
+ nBodyAssembler_.assembleTransferOperator();
+ nBodyAssembler_.assembleObstacle();
+
+ for (size_t i=0; i<nCouplings(); i++) {
+ auto& coupling = couplings_[i];
+ const auto& contactCoupling = nBodyAssembler_.getContactCouplings()[i]; // dual mortar object holding boundary patches
+ const auto nonmortarIdx = coupling->gridIdx_[0];
+ //const auto mortarIdx = coupling->gridIdx_[1];
+
+ frictionBoundaries_[nonmortarIdx]->addPatch(contactCoupling->nonmortarBoundary());
+ //frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
+ }
+
+ // assemble state energy norm
+ frictionBoundaryMass_.resize(nBodies());
+ for (size_t i=0; i<nBodies(); i++) {
+ frictionBoundaryMass_[i] = std::make_unique<ScalarMatrix>();
+ bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], *frictionBoundaryMass_[i]);
+ *frictionBoundaryMass_[i] /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
+
+ stateEnergyNorms_[i] = std::make_unique<typename LeafBody::StateEnergyNorm>(*frictionBoundaryMass_[i]);
+ }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::assembleFriction(
+ const Config::FrictionModel& frictionModel,
+ const std::vector<ScalarVector>& weightedNormalStress) {
+
+ assert(weightedNormalStress.size() == bodies_.size());
+
+ const size_t nBodies = bodies_.size();
+
+ // Lumping of the nonlinearity
+ std::vector<ScalarVector> weights(nBodies);
+ {
+ NeumannBoundaryAssembler<GridType, typename ScalarVector::block_type>
+ frictionalBoundaryAssembler(std::make_shared<
+ ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+ 1));
+
+ for (size_t i=0; i<nBodies; i++) {
+ bodies_[i]->assembler()->assembleBoundaryFunctional(frictionalBoundaryAssembler, weights[i], *frictionBoundaries_[i]);
+ }
+ }
+
+ /* globalFriction_ = std::make_shared<GlobalRateStateFriction<
+ Matrix, Vector, ZeroFunction, DeformedGridType>>(
+ nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress); */
+
+ switch (frictionModel) {
+ case Config::Truncated:
+ globalFriction_ = std::make_shared<GlobalRateStateFriction<
+ Matrix, VectorType, TruncatedRateState, GridType>>(
+ nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+ break;
+ case Config::Regularised:
+ globalFriction_ = std::make_shared<GlobalRateStateFriction<
+ Matrix, VectorType, RegularisedRateState, GridType>>(
+ nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+ break;
+ default:
+ assert(false);
+ break;
+ }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setBodies(const std::vector<std::shared_ptr<LeafBody>> 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++) {
+ matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
+ matrices_.damping[i] = bodies_[i]->matrices().damping;
+ matrices_.mass[i] = bodies_[i]->matrices().mass;
+ }
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
+ assert(this->nCouplings()==couplings.size());
+ couplings_ = couplings;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::constructBody(
+ const std::shared_ptr<BodyData<dim>>& bodyData,
+ const std::shared_ptr<HostGridType>& grid,
+ std::shared_ptr<LeafBody>& body) const {
+
+ body = std::make_shared<LeafBody>(bodyData, grid);
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::setDeformation(const std::vector<VectorType>& totalDeformation) {
+ assert(totalDeformation.size() == nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ body(i)->setDeformation(totalDeformation[i]);
+ }
+
+ nBodyAssembler_.assembleTransferOperator();
+ nBodyAssembler_.assembleObstacle();
+
for (size_t i=0; i<levelContactNetworks_.size(); i++) {
- levelContactNetworks_[i]->assemble();
+ levelContactNetworks_[i]->build();
}
}
-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 HostGridType, class VectorType>
+template <class LevelBoundaryPatch>
+void ContactNetwork<HostGridType, VectorType>::constructCoupling(
+ int nonMortarBodyIdx, int mortarBodyIdx,
+ const std::shared_ptr<LevelBoundaryPatch>& nonmortarPatch,
+ const std::shared_ptr<LevelBoundaryPatch>& mortarPatch,
+ const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+ const std::shared_ptr<GlobalFrictionData<dim>>& frictionData,
+ std::shared_ptr<FrictionCouplingPair>& coupling) const {
- levelContactNetworks_[level] = levelContactNetwork;
+ coupling = std::make_shared<FrictionCouplingPair>();
+
+ auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<BoundaryPatch>>();
+ std::shared_ptr<typename FrictionCouplingPair::BackEndType> backend = nullptr;
+ coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
+ coupling->setWeakeningPatch(weakeningPatch);
+ coupling->setFrictionData(frictionData);
}
-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());
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::level(size_t level) const
+-> const std::shared_ptr<LevelContactNetwork>& {
- levelContactNetworks_[level] = std::make_shared<LevelContactNetwork>(nBodies, nCouplings, level);
return levelContactNetworks_[level];
}
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nLevels() const
+-> size_t {
+
+ return levelContactNetworks_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodies() const
+-> size_t {
+
+ return bodies_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nCouplings() const
+-> size_t {
+
+ return couplings_.size();
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::body(int i) const
+-> const std::shared_ptr<LeafBody>& {
+
+ return bodies_[i];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::coupling(int i) const
+-> const std::shared_ptr<FrictionCouplingPair>& {
+
+ return couplings_[i];
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::couplings() const
+-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
+
+ return couplings_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler()
+-> Dune::Contact::NBodyAssembler<GridType, VectorType>& {
+
+ return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::nBodyAssembler() const
+-> const Dune::Contact::NBodyAssembler<GridType, VectorType>& {
+
+ return nBodyAssembler_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::matrices() const
+-> const Matrices<Matrix,2>& {
+
+ return matrices_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::stateEnergyNorms() const
+-> const StateEnergyNorms& {
+
+ return stateEnergyNorms_;
+}
+
+template <class HostGridType, class VectorType>
+auto ContactNetwork<HostGridType, VectorType>::globalFriction() const
+-> GlobalFriction& {
+
+ return *globalFriction_;
+}
+
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::externalForces(ExternalForces& externalForces) const {
+ externalForces.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ externalForces[i] = std::make_unique<typename LeafBody::ExternalForce>(bodies_[i]->externalForce());
+ }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::totalNodes(
+ const std::string& tag,
+ Dune::BitSetVector<dim>& totalNodes) const {
+
+ int totalSize = 0;
+ for (size_t i=0; i<nBodies(); i++) {
+ totalSize += this->body(i)->nVertices();
+ }
+
+ 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 LeafBody::BoundaryCondition>> boundaryConditions;
+ body->boundaryConditions(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<dim; k++)
+ totalNodes[idx][k] = (*nodes)[j][k];
+ idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+ }
+ }
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a vector of BitSetVector
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::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 HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::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 HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::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 HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::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 HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const {
+ nodes.resize(nBodies());
+ for (size_t i=0; i<nBodies(); i++) {
+ nodes[i] = frictionBoundaries_[i]->getVertices();
+ }
+}
+
#include "contactnetwork_tmpl.cc"
diff --git a/src/data-structures/contactnetwork.hh b/src/data-structures/contactnetwork.hh
index 25baf45b..6fb2e48b 100644
--- a/src/data-structures/contactnetwork.hh
+++ b/src/data-structures/contactnetwork.hh
@@ -1,37 +1,146 @@
#ifndef SRC_CONTACTNETWORK_HH
#define SRC_CONTACTNETWORK_HH
+#include <memory>
+
+#include <dune/common/promotiontraits.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+#include <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "../assemblers.hh"
+#include "../frictioncouplingpair.hh"
+#include "body.hh"
+#include "enums.hh"
+#include "matrices.hh"
#include "levelcontactnetwork.hh"
-template <class GridType, int dim> class ContactNetwork {
- public:
- using LevelContactNetwork = LevelContactNetwork<GridType, dim>;
+template <class HostGridType, class VectorType>
+class ContactNetwork {
+ //using Vector = Dune::BlockVector<Dune::FieldVector<double, dim>>;
+public:
+ enum {dim = HostGridType::dimensionworld};
+
+ using field_type = typename Dune::PromotionTraits<typename VectorType::field_type,
+ typename HostGridType::ctype>::PromotedType;
+
+ using LeafBody = LeafBody<HostGridType, VectorType>;
+ using GridType = typename LeafBody::GridType;
+ using GridView = typename LeafBody::GridView;
+ using BoundaryPatch = BoundaryPatch<GridView>;
+
+ using Assembler = typename LeafBody::Assembler;
+ using Matrix = typename Assembler::Matrix;
+
+ using ScalarMatrix = typename Assembler::ScalarMatrix;
+ using ScalarVector = typename Assembler::ScalarVector;
+
+ using LocalVector = typename VectorType::block_type;
+ using FrictionCouplingPair = FrictionCouplingPair<GridType, LocalVector, field_type>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+
+ using BoundaryFunctions = std::vector<typename LeafBody::BoundaryFunctions>;
+ using BoundaryNodes = std::vector<typename LeafBody::BoundaryNodes>;
+ using BoundaryPatchNodes = std::vector<typename LeafBody::BoundaryPatchNodes>;
+ using BoundaryPatches = std::vector<typename LeafBody::BoundaryPatches>;
+
+ using GlobalFriction = GlobalFriction<Matrix, VectorType>;
+
+ using StateEnergyNorms = std::vector<std::unique_ptr<typename LeafBody::StateEnergyNorm> >;
+ using ExternalForces = std::vector<std::unique_ptr<const typename LeafBody::ExternalForce> >;
+
+ using NBodyAssembler = Dune::Contact::NBodyAssembler<GridType, VectorType>;
+
+ using LevelContactNetwork = LevelContactNetwork<GridType, FrictionCouplingPair, field_type>;
+ using Body = typename LevelContactNetwork::Body;
+
+private:
+ std::vector<std::shared_ptr<LevelContactNetwork>> levelContactNetworks_;
+
+ std::vector<std::shared_ptr<LeafBody>> bodies_;
+ std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+ std::vector<std::unique_ptr<BoundaryPatch>> frictionBoundaries_; // union of all boundary patches per body
+ std::vector<std::unique_ptr<ScalarMatrix>> frictionBoundaryMass_;
+ StateEnergyNorms stateEnergyNorms_;
+
+ NBodyAssembler nBodyAssembler_;
+
+ Matrices<Matrix,2> matrices_;
+
+ std::shared_ptr<GlobalFriction> globalFriction_;
+
+public:
+ ContactNetwork(size_t nBodies, size_t nCouplings, size_t nLevels = 0);
+
+ void addLevel(std::shared_ptr<LevelContactNetwork> levelContactNetwork);
+ void addLevel(const std::vector<size_t>& bodyLevels, const size_t level);
+
+ void assemble();
+ void assembleFriction(
+ const Config::FrictionModel& frictionModel,
+ const std::vector<ScalarVector>& weightedNormalStress);
+
+ void setBodies(const std::vector<std::shared_ptr<LeafBody>> bodies);
+ void setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings);
+
+ void setDeformation(const std::vector<VectorType>& totalDeformation);
+
+ void constructBody(
+ const std::shared_ptr<BodyData<dim>>& bodyData,
+ const std::shared_ptr<HostGridType>& grid,
+ std::shared_ptr<LeafBody>& body) const;
+
+ template <class LevelBoundaryPatch>
+ void constructCoupling(
+ int nonMortarBodyIdx,
+ int mortarBodyIdx,
+ const std::shared_ptr<LevelBoundaryPatch>& nonmortarPatch,
+ const std::shared_ptr<LevelBoundaryPatch>& mortarPatch,
+ const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
+ const std::shared_ptr<GlobalFrictionData<dim>>& frictionData,
+ std::shared_ptr<FrictionCouplingPair>& coupling) const;
- private:
- std::vector<std::shared_ptr<LevelContactNetwork>> levelContactNetworks_;
+ // getter
+ auto level(size_t level) const -> const std::shared_ptr<LevelContactNetwork>&;
- public:
- ContactNetwork() {}
+ auto nLevels() const -> size_t;
+ auto nBodies() const -> size_t;
+ auto nCouplings() const -> size_t;
- void addLevelContactNetwork(std::shared_ptr<LevelContactNetwork> levelContactNetwork);
- const std::shared_ptr<LevelContactNetwork>& addLevelContactNetwork(int nBodies, int nCouplings, int level);
+ auto body(int i) const -> const std::shared_ptr<LeafBody>&;
- void assemble();
+ auto coupling(int i) const -> const std::shared_ptr<FrictionCouplingPair>&;
+ auto couplings() const -> const std::vector<std::shared_ptr<FrictionCouplingPair>>&;
- const std::shared_ptr<LevelContactNetwork>& levelContactNetwork(size_t level) const {
- return levelContactNetworks_[level];
- }
+ auto nBodyAssembler() -> NBodyAssembler&;
+ auto nBodyAssembler() const -> const NBodyAssembler&;
- void resize(size_t size) {
- levelContactNetworks_.resize(size);
- }
+ auto matrices() const -> const Matrices<Matrix,2>&;
- size_t size() const {
- return levelContactNetworks_.size();
- }
+ auto stateEnergyNorms() const -> const StateEnergyNorms&;
+ auto globalFriction() const -> GlobalFriction&;
+ void externalForces(ExternalForces& externalForces) const;
- size_t maxLevel() const {
- return size()-1;
- }
+ void totalNodes(
+ const std::string& tag,
+ Dune::BitSetVector<dim>& 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 frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const;
};
#endif
diff --git a/src/data-structures/contactnetwork_tmpl.cc b/src/data-structures/contactnetwork_tmpl.cc
index 32c9a280..03374e69 100644
--- a/src/data-structures/contactnetwork_tmpl.cc
+++ b/src/data-structures/contactnetwork_tmpl.cc
@@ -3,5 +3,10 @@
#endif
#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
-template class ContactNetwork<Grid, MY_DIM>;
+#include "contactnetwork.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+template class ContactNetwork<Grid, Vector>;
diff --git a/src/data-structures/levelcontactnetwork.cc b/src/data-structures/levelcontactnetwork.cc
index 5edb70fc..f7a34d4a 100644
--- a/src/data-structures/levelcontactnetwork.cc
+++ b/src/data-structures/levelcontactnetwork.cc
@@ -6,6 +6,7 @@
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/facehierarchy.hh>
#include <dune/tectonic/globalratestatefriction.hh>
#include <dune/tectonic/frictionpotential.hh>
@@ -16,186 +17,51 @@
#include "../utils/tobool.hh"
#include "../utils/debugutils.hh"
-template <class GridType, int dimension>
-LevelContactNetwork<GridType, dimension>::LevelContactNetwork(
+template <class GridType, class FrictionCouplingPair, class field_type>
+LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::LevelContactNetwork(
int nBodies,
int nCouplings,
int level) :
level_(level),
bodies_(nBodies),
couplings_(nCouplings),
- deformedGrids_(nBodies),
- leafViews_(nBodies),
- levelViews_(nBodies),
- leafVertexCounts_(nBodies),
- frictionBoundaries_(nCouplings),
- stateEnergyNorms_(nBodies),
- nBodyAssembler_(nBodies, nCouplings)
+ nonmortarPatches_(nCouplings),
+ mortarPatches_(nCouplings),
+ glues_(nCouplings)
{}
-template <class GridType, int dimension>
-LevelContactNetwork<GridType, dimension>::~LevelContactNetwork() {
- for (size_t i=0; i<nBodies(); i++) {
- delete frictionBoundaries_[i];
- delete frictionBoundaryMass_[i];
- delete stateEnergyNorms_[i];
- }
-}
-
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::assemble() {
- using BoundaryPatch = DeformedLeafBoundaryPatch;
-
- for (size_t i=0; i<nBodies(); i++) {
- bodies_[i]->assemble();
- //printBasisDofLocation(bodies_[i]->assembler()->vertexBasis); //TODO remove after debugging
-
- frictionBoundaries_[i] = new BoundaryPatch(bodies_[i]->leafView(), false);
- }
-
- // set up dune-contact nBodyAssembler
- nBodyAssembler_.setGrids(deformedGrids_);
-
- for (size_t i=0; i<nCouplings(); i++) {
- nBodyAssembler_.setCoupling(*couplings_[i], i);
- }
-
- nBodyAssembler_.assembleTransferOperator();
- nBodyAssembler_.assembleObstacle();
-
- for (size_t i=0; i<nCouplings(); i++) {
- auto& coupling = couplings_[i];
- const auto& contactCoupling = nBodyAssembler_.getContactCouplings()[i]; // dual mortar object holding boundary patches
- const auto nonmortarIdx = coupling->gridIdx_[0];
- //const auto mortarIdx = coupling->gridIdx_[1];
-
- frictionBoundaries_[nonmortarIdx]->addPatch(contactCoupling->nonmortarBoundary());
- //frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
- }
-
- // assemble state energy norm
- frictionBoundaryMass_.resize(nBodies());
- for (size_t i=0; i<nBodies(); i++) {
- frictionBoundaryMass_[i] = new ScalarMatrix();
- bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], *frictionBoundaryMass_[i]);
- *frictionBoundaryMass_[i] /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
- stateEnergyNorms_[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(*frictionBoundaryMass_[i]);
- }
-}
-
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::assembleFriction(
- const Config::FrictionModel& frictionModel,
- const std::vector<ScalarVector>& weightedNormalStress) {
-
- assert(weightedNormalStress.size() == bodies_.size());
-
- const size_t nBodies = bodies_.size();
-
- // Lumping of the nonlinearity
- std::vector<ScalarVector> weights(nBodies);
- {
- NeumannBoundaryAssembler<DeformedGridType, typename ScalarVector::block_type>
- frictionalBoundaryAssembler(std::make_shared<
- ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
- 1));
-
- for (size_t i=0; i<nBodies; i++) {
- bodies_[i]->assembler()->assembleBoundaryFunctional(frictionalBoundaryAssembler, weights[i], *frictionBoundaries_[i]);
- }
- }
-
- /* globalFriction_ = std::make_shared<GlobalRateStateFriction<
- Matrix, Vector, ZeroFunction, DeformedGridType>>(
- nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress); */
-
- switch (frictionModel) {
- case Config::Truncated:
- globalFriction_ = std::make_shared<GlobalRateStateFriction<
- Matrix, Vector, TruncatedRateState, DeformedGridType>>(
- nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
- break;
- case Config::Regularised:
- globalFriction_ = std::make_shared<GlobalRateStateFriction<
- Matrix, Vector, RegularisedRateState, DeformedGridType>>(
- nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
- break;
- default:
- assert(false);
- break;
- }
-}
-
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::shared_ptr<Body>> bodies) {
- assert(nBodies()==bodies.size());
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::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].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);
-
- 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 dimension>
-void LevelContactNetwork<GridType, dimension>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
- assert(this->nCouplings()==couplings.size());
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings) {
+ assert(nCouplings() == couplings.size());
couplings_ = couplings;
}
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::constructBody(
- const std::shared_ptr<BodyData<dimension>>& bodyData,
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::constructBody(
+ const std::shared_ptr<BodyData<dimworld>>& bodyData,
const std::shared_ptr<GridType>& grid,
+ const size_t level,
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,
- const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
- const std::shared_ptr<GlobalFrictionData<dimension>>& frictionData,
- std::shared_ptr<FrictionCouplingPair>& coupling) const {
-
- coupling = std::make_shared<FrictionCouplingPair>();
-
- auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<DeformedLeafBoundaryPatch>>();
- std::shared_ptr<typename FrictionCouplingPair::BackEndType> backend = nullptr;
- coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
- coupling->setWeakeningPatch(weakeningPatch);
- coupling->setFrictionData(frictionData);
+ body = std::make_shared<Body>(bodyData, grid, level);
}
// 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(
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::totalNodes(
const std::string& tag,
- Dune::BitSetVector<dimension>& totalNodes) const {
+ Dune::BitSetVector<dimworld>& totalNodes) const {
int totalSize = 0;
for (size_t i=0; i<nBodies(); i++) {
- totalSize += this->body(i)->leafVertexCount();
+ totalSize += this->body(i)->nVertices();
}
totalNodes.resize(totalSize);
@@ -203,14 +69,14 @@ void LevelContactNetwork<GridType, dimension>::totalNodes(
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);
+ std::vector<std::shared_ptr<typename Body::BoundaryCondition>> boundaryConditions;
+ body->boundaryConditions(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++)
+ for (int k=0; k<dimworld; k++)
totalNodes[idx][k] = (*nodes)[j][k];
idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
}
@@ -218,8 +84,8 @@ void LevelContactNetwork<GridType, dimension>::totalNodes(
}
// 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(
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryPatchNodes(
const std::string& tag,
BoundaryPatchNodes& nodes) const {
@@ -231,8 +97,8 @@ void LevelContactNetwork<GridType, dimension>::boundaryPatchNodes(
}
// 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(
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryNodes(
const std::string& tag,
BoundaryNodes& nodes) const {
@@ -244,8 +110,8 @@ void LevelContactNetwork<GridType, dimension>::boundaryNodes(
}
// collects all leaf boundary functions from the different bodies identified by "tag"
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::boundaryFunctions(
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryFunctions(
const std::string& tag,
BoundaryFunctions& functions) const {
@@ -257,8 +123,8 @@ void LevelContactNetwork<GridType, dimension>::boundaryFunctions(
}
// collects all leaf boundary patches from the different bodies identified by "tag"
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::boundaryPatches(
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::boundaryPatches(
const std::string& tag,
BoundaryPatches& patches) const {
@@ -269,141 +135,113 @@ void LevelContactNetwork<GridType, dimension>::boundaryPatches(
}
}
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const {
- nodes.resize(nBodies());
- for (size_t i=0; i<nBodies(); i++) {
- nodes[i] = frictionBoundaries_[i]->getVertices();
- }
-}
-
-
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::stateEnergyNorms() const
--> const StateEnergyNorms& {
-
- return stateEnergyNorms_;
-}
-
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::level() const
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::level() const
-> int {
return level_;
}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::nBodies() const
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::nBodies() const
-> size_t {
return bodies_.size();
}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::nCouplings() const
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::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
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::body(int i) const
-> const std::shared_ptr<Body>& {
return bodies_[i];
}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::coupling(int i) const
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::coupling(int i) const
-> const std::shared_ptr<FrictionCouplingPair>& {
return couplings_[i];
}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::deformedGrids() const
--> const std::vector<const DeformedGridType*>& {
-
- return deformedGrids_;
-}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::couplings() const
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::couplings() const
-> const std::vector<std::shared_ptr<FrictionCouplingPair>>& {
return couplings_;
}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::leafView(int bodyID) const
--> const DeformedLeafGridView& {
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::prolong(const BoundaryPatch& coarsePatch, BoundaryPatch& finePatch, const size_t fineLevel) {
- return *leafViews_[bodyID];
-}
+ if (not(coarsePatch.isInitialized()))
+ DUNE_THROW(Dune::Exception, "Coarse boundary patch has not been set up correctly!");
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::levelViews(int bodyID) const
--> const std::vector<std::unique_ptr<const DeformedLevelGridView>>& {
+ const GridType& grid = coarsePatch.gridView().grid();
- return levelViews_[bodyID];
-}
+ // Set array sizes correctly
+ finePatch.setup(grid.levelGridView(fineLevel));
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::levelView(int bodyID, int level) const
--> const DeformedLevelGridView& {
+ for (const auto& pIt : coarsePatch) {
+ const auto& inside = pIt.inside();
- return *levelViews_[bodyID][level];
-}
+ if (inside.level() == fineLevel)
+ finePatch.addFace(inside, pIt.indexInInside());
+ else {
+ Face<GridType> face(grid, inside, pIt.indexInInside());
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::leafVertexCount(int bodyID) const
--> int {
+ typename Face<GridType>::HierarchicIterator it = face.hbegin(fineLevel);
+ typename Face<GridType>::HierarchicIterator end = face.hend(fineLevel);
- return leafVertexCounts_[bodyID];
+ for(; it!=end; ++it) {
+ if (it->element().level() == fineLevel)
+ finePatch.addFace(it->element(), it->index());
+ }
+ }
+ }
}
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::nBodyAssembler()
--> Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& {
-
- return nBodyAssembler_;
-}
+template <class GridType, class FrictionCouplingPair, class field_type>
+void LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::build(field_type overlap) {
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::nBodyAssembler() const
--> const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& {
+ for (size_t i=0; i<couplings_.size(); i++) {
+ auto& coupling = *couplings_[i];
- return nBodyAssembler_;
-}
+ const size_t nmBody = coupling.gridIdx_[0];
+ const size_t mBody = coupling.gridIdx_[1];
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::matrices() const
--> const Matrices<Matrix,2>& {
+ this->prolong(*coupling.patch0(), nonmortarPatches_[i], body(nmBody)->level());
+ this->prolong(*coupling.patch1(), mortarPatches_[i], body(mBody)->level());
- return matrices_;
-}
+ using Element = typename GridView::template Codim<0>::Entity;
+ auto desc0 = [&] (const Element& e, unsigned int face) {
+ return nonmortarPatches_[i].contains(e,face);
+ };
+ auto desc1 = [&] (const Element& e, unsigned int face) {
+ return mortarPatches_[i].contains(e,face);
+ };
-template <class GridType, int dimension>
-auto LevelContactNetwork<GridType, dimension>::globalFriction() const
--> GlobalFriction& {
+ auto extract0 = std::make_shared<Extractor>(body(nmBody)->gridView(), desc0);
+ auto extract1 = std::make_shared<Extractor>(body(mBody)->gridView(), desc1);
- return *globalFriction_;
+ gridGlueBackend_ = std::make_shared< Dune::GridGlue::ContactMerge<dim, ctype> >(overlap);
+ glues_[i] = std::make_shared<Glue>(extract0, extract1, gridGlueBackend_);
+ glues_[i]->build();
+ }
}
-template <class GridType, int dimension>
-void LevelContactNetwork<GridType, dimension>::externalForces(ExternalForces& externalForces) const {
- externalForces.resize(nBodies());
+template <class GridType, class FrictionCouplingPair, class field_type>
+auto LevelContactNetwork<GridType, FrictionCouplingPair, field_type>::glue(int i) const
+-> const std::shared_ptr<Glue>& {
- for (size_t i=0; i<nBodies(); i++) {
- externalForces[i] = &bodies_[i]->externalForce();
- }
+ return glues_[i];
}
#include "levelcontactnetwork_tmpl.cc"
diff --git a/src/data-structures/levelcontactnetwork.hh b/src/data-structures/levelcontactnetwork.hh
index 68bb716a..c4976718 100644
--- a/src/data-structures/levelcontactnetwork.hh
+++ b/src/data-structures/levelcontactnetwork.hh
@@ -6,12 +6,18 @@
#include <dune/istl/bvector.hh>
+#include <dune/grid-glue/gridglue.hh>
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+#include <dune/grid-glue/merging/merger.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/fufem/boundarypatch.hh>
+
#include <dune/tectonic/bodydata.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/globalfrictiondata.hh>
@@ -23,134 +29,84 @@
#include "matrices.hh"
//#include "multi-body-problem-data/myglobalfrictiondata.hh"
-template <class GridType, int dimension> class LevelContactNetwork {
- private:
- //static const int dimension = typename GridType::dimension;
- using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
-
- 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 FrictionCouplingPair = FrictionCouplingPair<DeformedGridType, LocalVector, 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 GlobalFriction = GlobalFriction<Matrix, Vector>;
-
- using StateEnergyNorms = std::vector<const typename Body::StateEnergyNorm *>;
- using ExternalForces = std::vector<const typename Body::ExternalForce *>;
-
- using NBodyAssembler = Dune::Contact::NBodyAssembler<DeformedGridType, Vector>;
-
- public:
- LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
- ~LevelContactNetwork();
-
- void assemble();
- void assembleFriction(
- const Config::FrictionModel& frictionModel,
- const std::vector<ScalarVector>& weightedNormalStress);
-
- void setBodies(const std::vector<std::shared_ptr<Body>> bodies);
- void setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> 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,
- const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch,
- const std::shared_ptr<GlobalFrictionData<dimension>>& frictionData,
- std::shared_ptr<FrictionCouplingPair>& 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 frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const;
-
- auto stateEnergyNorms() const -> const StateEnergyNorms&;
- auto globalFriction() const -> GlobalFriction&;
- void externalForces(ExternalForces& externalForces) const;
-
- auto level() const -> int;
-
- auto nBodies() const -> size_t;
- auto nCouplings() const -> size_t;
-
- auto body(int i) const -> const std::shared_ptr<Body>&;
-
- auto coupling(int i) const -> const std::shared_ptr<FrictionCouplingPair>&;
- auto couplings() const -> const std::vector<std::shared_ptr<FrictionCouplingPair>>&;
-
- auto deformedGrids() const -> const std::vector<const DeformedGridType*>& ;
-
- auto leafVertexCounts() const -> const std::vector<size_t>&;
- auto leafVertexCount(int bodyID) const -> int;
-
- 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>>&;
-
- auto nBodyAssembler() -> Dune::Contact::NBodyAssembler<DeformedGridType, Vector>&;
- auto nBodyAssembler() const -> const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>&;
-
- auto matrices() const -> const Matrices<Matrix,2>&;
-
- private:
- const int level_;
-
- std::vector<std::shared_ptr<Body>> bodies_;
- std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
-
- std::vector<const DeformedGridType*> deformedGrids_;
- std::vector<std::unique_ptr<const DeformedLeafGridView>> leafViews_;
- std::vector<std::vector<std::unique_ptr<const DeformedLevelGridView>>> levelViews_;
- std::vector<size_t> leafVertexCounts_;
-
- std::vector<DeformedLeafBoundaryPatch*> frictionBoundaries_; // union of all boundary patches per body
- std::vector<ScalarMatrix*> frictionBoundaryMass_;
- StateEnergyNorms stateEnergyNorms_;
-
- NBodyAssembler nBodyAssembler_;
-
- Matrices<Matrix,2> matrices_;
-
- std::shared_ptr<GlobalFriction> globalFriction_;
+template <class GridTypeTEMPLATE, class FrictionCouplingPair, class field_type>
+class LevelContactNetwork {
+public:
+ using GridType = GridTypeTEMPLATE;
+
+ enum {dim = GridType::dimension};
+ enum {dimworld = GridType::dimensionworld};
+
+ using ctype = typename GridType::ctype;
+
+ using GridView = typename GridType::LevelGridView;
+ using Body = Body<GridView>;
+
+ 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 BoundaryPatch = BoundaryPatch<GridView>;
+
+ using Extractor = Dune::GridGlue::Codim1Extractor<GridView>;
+ using Glue = Dune::GridGlue::GridGlue<Extractor, Extractor>;
+
+public:
+ LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
+
+ void setBodies(const std::vector<std::shared_ptr<Body>> bodies);
+ void setCouplings(const std::vector<std::shared_ptr<FrictionCouplingPair>> couplings);
+
+ void constructBody(
+ const std::shared_ptr<BodyData<dimworld>>& bodyData,
+ const std::shared_ptr<GridType>& grid,
+ const size_t level,
+ std::shared_ptr<Body>& body) const;
+
+ void build(field_type overlap = 1e-2);
+
+ // getter
+ void totalNodes(
+ const std::string& tag,
+ Dune::BitSetVector<dimworld>& 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;
+
+ auto level() const -> int;
+
+ auto nBodies() const -> size_t;
+ auto nCouplings() const -> size_t;
+
+ auto body(int i) const -> const std::shared_ptr<Body>&;
+
+ auto coupling(int i) const -> const std::shared_ptr<FrictionCouplingPair>&;
+ auto couplings() const -> const std::vector<std::shared_ptr<FrictionCouplingPair>>&;
+
+ auto glue(int i) const -> const std::shared_ptr<Glue>&;
+
+private:
+ void prolong(const BoundaryPatch& coarsePatch, BoundaryPatch& finePatch, const size_t fineLevel);
+
+ const int level_;
+
+ std::vector<std::shared_ptr<Body>> bodies_;
+ std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+ std::vector<BoundaryPatch> nonmortarPatches_;
+ std::vector<BoundaryPatch> mortarPatches_;
+ std::shared_ptr< Dune::GridGlue::Merger<field_type, dim-1, dim-1, dim> > gridGlueBackend_;
+ std::vector<std::shared_ptr<Glue>> glues_;
};
#endif
diff --git a/src/data-structures/levelcontactnetwork_tmpl.cc b/src/data-structures/levelcontactnetwork_tmpl.cc
index 1bd54512..da253a7f 100644
--- a/src/data-structures/levelcontactnetwork_tmpl.cc
+++ b/src/data-structures/levelcontactnetwork_tmpl.cc
@@ -3,8 +3,10 @@
#endif
#include "../explicitgrid.hh"
-#include "levelcontactnetwork.hh"
+#include "../explicitvectors.hh"
-using MyLevelContactNetwork = LevelContactNetwork<Grid, MY_DIM>;
+#include "../frictioncouplingpair.hh"
+#include "contactnetwork_tmpl.cc"
+#include "levelcontactnetwork.hh"
-template class LevelContactNetwork<Grid, MY_DIM>;
+template class LevelContactNetwork<typename MyContactNetwork::GridType, typename MyContactNetwork::FrictionCouplingPair, typename MyContactNetwork::field_type>;
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index bab650b9..9c5dd4af 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -13,7 +13,7 @@
#include <dune/tectonic/bodydata.hh>
#include "../assemblers.hh"
-#include "levelcontactnetwork.hh"
+#include "contactnetwork.hh"
#include "matrices.hh"
#include "../spatial-solving/solverfactory.hh"
#include "../spatial-solving/tnnmg/functional.hh"
@@ -86,9 +86,9 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Set up initial conditions
template <class GridType>
- void setupInitialConditions(const Dune::ParameterTree& parset, const LevelContactNetwork<GridType, dims>& levelContactNetwork) {
- using Matrix = typename LevelContactNetwork<GridType, dims>::Matrix;
- const auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
+ void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork<GridType, Vector>& contactNetwork) {
+ using Matrix = typename ContactNetwork<GridType, Vector>::Matrix;
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
@@ -227,14 +227,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]);
+ contactNetwork.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.totalNodes("dirichlet", dirichletNodes);
+ contactNetwork.totalNodes("dirichlet", dirichletNodes);
/*for (size_t i=0; i<dirichletNodes.size(); i++) {
bool val = false;
for (size_t d=0; d<dims; d++) {
@@ -247,7 +247,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
}
}*/
- solveLinearProblem(dirichletNodes, levelContactNetwork.matrices().elasticity, ell0, u,
+ solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u,
parset.sub("u0.solver"));
print(u, "initial u:");
@@ -261,9 +261,9 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Initial normal stress
- const auto& body = levelContactNetwork.body(i);
- std::vector<std::shared_ptr<typename LevelContactNetwork<GridType, dims>::Body::LeafBoundaryCondition>> frictionBoundaryConditions;
- body->leafBoundaryConditions("friction", frictionBoundaryConditions);
+ const auto& body = contactNetwork.body(i);
+ std::vector<std::shared_ptr<typename ContactNetwork<GridType, Vector>::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
+ body->boundaryConditions("friction", frictionBoundaryConditions);
for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
@@ -278,7 +278,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
}
Dune::BitSetVector<dims> noNodes(dirichletNodes.size(), false);
- solveLinearProblem(noNodes, levelContactNetwork.matrices().mass, accelerationRHS, a,
+ solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
parset.sub("a0.solver"));
print(a, "initial a:");
diff --git a/src/explicitvectors.hh b/src/explicitvectors.hh
index 548c681f..54887ada 100644
--- a/src/explicitvectors.hh
+++ b/src/explicitvectors.hh
@@ -7,12 +7,14 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
-using LocalVector = Dune::FieldVector<double, MY_DIM>;
-using LocalMatrix = Dune::FieldMatrix<double, MY_DIM, MY_DIM>;
+using ctype = double;
+
+using LocalVector = Dune::FieldVector<ctype, MY_DIM>;
+using LocalMatrix = Dune::FieldMatrix<ctype, MY_DIM, MY_DIM>;
using Vector = Dune::BlockVector<LocalVector>;
using Matrix = Dune::BCRSMatrix<LocalMatrix>;
-using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+using ScalarVector = Dune::BlockVector<Dune::FieldVector<ctype, 1>>;
+using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<ctype, 1, 1>>;
using BitVector = Dune::BitSetVector<MY_DIM>;
#endif
diff --git a/src/factories/contactnetworkfactory.hh b/src/factories/contactnetworkfactory.hh
index 90298ad6..c1968822 100644
--- a/src/factories/contactnetworkfactory.hh
+++ b/src/factories/contactnetworkfactory.hh
@@ -4,45 +4,55 @@
#include <dune/common/parametertree.hh>
#include "../data-structures/contactnetwork.hh"
-#include "../data-structures/levelcontactnetwork.hh"
-#include "levelcontactnetworkfactory.hh"
-
-template <class GridType, int dim>
+template <class HostGridType, class VectorType>
class ContactNetworkFactory {
public:
- using ContactNetwork = ContactNetwork<GridType, dim>;
+ using ContactNetwork = ContactNetwork<HostGridType, VectorType>;
protected:
- using LevelContactNetworkFactory = LevelContactNetworkFactory<GridType, dim>;
+ using LeafBody = typename ContactNetwork::LeafBody;
+ using FrictionCouplingPair = typename ContactNetwork::FrictionCouplingPair;
const Dune::ParameterTree& parset_;
+ const size_t bodyCount_;
+ const size_t couplingCount_;
- std::vector<std::shared_ptr<LevelContactNetworkFactory>> levelContactNetworkFactories_;
+ std::vector<std::shared_ptr<LeafBody>> bodies_;
+ std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
ContactNetwork contactNetwork_;
private:
virtual void setBodies() = 0;
+ virtual void setLevelBodies() = 0;
virtual void setCouplings() = 0;
+ virtual void setLevelCouplings() = 0;
virtual void setBoundaryConditions() = 0;
public:
- ContactNetworkFactory(const Dune::ParameterTree& parset) :
+ ContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
parset_(parset),
- contactNetwork_() {}
+ bodyCount_(bodyCount),
+ couplingCount_(couplingCount),
+ bodies_(bodyCount_),
+ couplings_(couplingCount_),
+ contactNetwork_(bodyCount_, couplingCount_) {}
void build() {
setBodies();
- contactNetwork_.setBodies();
+ contactNetwork_.setBodies(bodies_);
setCouplings();
- contactNetwork_.setCouplings();
+ contactNetwork_.setCouplings(couplings_);
+
+ setLevelBodies();
+ setLevelCouplings();
- setBoundaryConditions();
+ setBoundaryConditions();
}
- ContactNetwork& levelContactNetwork() {
+ ContactNetwork& contactNetwork() {
return contactNetwork_;
}
};
diff --git a/src/factories/levelcontactnetworkfactory.hh b/src/factories/levelcontactnetworkfactory.hh
index 833036c6..d85b0853 100644
--- a/src/factories/levelcontactnetworkfactory.hh
+++ b/src/factories/levelcontactnetworkfactory.hh
@@ -1,12 +1,12 @@
-#ifndef SRC_LEVELCONTACTNETWORKFACTORY_HH
-#define SRC_LEVELCONTACTNETWORKFACTORY_HH
+#ifndef SRC_CONTACTNETWORKFACTORY_HH
+#define SRC_CONTACTNETWORKFACTORY_HH
#include <dune/common/parametertree.hh>
#include "../data-structures/levelcontactnetwork.hh"
-template <class GridType, int dims>
-class LevelContactNetworkFactory {
+template <class HostGridType>
+class ContactNetworkFactory {
public:
using LevelContactNetwork = LevelContactNetwork<GridType, dims>;
diff --git a/src/factories/stackedblocksfactory.cc b/src/factories/stackedblocksfactory.cc
index feb057cd..5a2aaaff 100644
--- a/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
@@ -9,7 +9,6 @@
#include "../multi-body-problem-data/bc.hh"
#include "../multi-body-problem-data/grid/cuboidgeometry.hh"
#include "../multi-body-problem-data/myglobalfrictiondata.hh"
-#include "../multi-body-problem-data/weakpatch.hh"
#include "../utils/diameter.hh"
diff --git a/src/factories/threeblocksfactory.cc b/src/factories/threeblocksfactory.cc
new file mode 100644
index 00000000..7dded09d
--- /dev/null
+++ b/src/factories/threeblocksfactory.cc
@@ -0,0 +1,204 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../multi-body-problem-data/bc.hh"
+#include "../multi-body-problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "threeblocksfactory.hh"
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
+ // set up cuboid geometries
+
+ std::array<double, 3> lengths = {3.0, 1.0, 1.0};
+ std::array<double, 3> heights = {1.0, 1.0, 1.0};
+
+ std::array<GlobalCoords, 3> origins;
+
+ const auto& subParset = this->parset_.sub("boundary.friction.weakening");
+
+#if MY_DIM == 3 // TODO: not implemented
+ double const depth = 0.60;
+
+ origins[0] = {0, 0, 0};
+ origins[1] = {0, origins[0][1] + widths[0], 0};
+ origins[2] = {origins[1][0] + lengths[1], origins[0][1] + widths[0], 0};
+
+ cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0], depth);
+ cuboidGeometries_[0]->addWeakeningPatch(subParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
+
+ cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1], depth);
+ cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
+ cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1] + heights[1], 0});
+
+ cuboidGeometries_[2] = std::make_shared<CuboidGeometry>(origins[2], lengths[2], heights[2], depth);
+ cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0] + lengths[2], origins[2][1], 0});
+ cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0], origins[2][1] + widths[2], 0});
+#elif MY_DIM == 2
+ origins[0] = {0, 0};
+ origins[1] = {0, origins[0][1] + heights[0]};
+ origins[2] = {origins[1][0] + lengths[1], origins[0][1] + heights[0]};
+
+ cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0]);
+ cuboidGeometries_[0]->addWeakeningPatch(subParset, {origins[0][0], origins[0][1]+ heights[0]}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0]});
+
+ cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1]);
+ cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1]});
+ cuboidGeometries_[1]->addWeakeningPatch(subParset, {origins[1][0] + lengths[1], origins[1][1]}, {origins[1][0] + lengths[1], origins[1][1] + heights[1]});
+
+ cuboidGeometries_[2] = std::make_shared<CuboidGeometry>(origins[2], lengths[2], heights[2]);
+ cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0] + lengths[2], origins[2][1]});
+ cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0], origins[2][1] + heights[2]});
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+ // set up reference grids
+ gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
+ auto& grids = gridConstructor_->getGrids();
+
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+ // define weak patch and refine grid
+ const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
+ for (size_t j=0; j<weakeningRegions.size(); j++) {
+ refine(*grids[i], weakeningRegions[j], 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++) {
+ this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, grids[i]);
+ }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
+ const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel(), this->bodies_[2]->grid()->maxLevel()});
+
+ for (size_t l=0; l<=maxLevel; l++) {
+ std::vector<size_t> bodyLevels(3, l);
+ this->contactNetwork_.addLevel(bodyLevels, l);
+ }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& cuboidGeometry = *cuboidGeometries_[i];
+ std::cout << this->bodies_[i]->nVertices() << std::endl;
+ std::cout << (leafFaces_[i] == nullptr) << std::endl;
+ leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
+ levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
+ }
+
+ auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+ std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
+
+ std::array<std::array<size_t, 2>, 3> couplingIndices;
+
+ nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
+ mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
+ weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[0]->weakeningRegions()[0]);
+ couplingIndices[0] = {0, 1};
+
+ nonmortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->lower);
+ mortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
+ weakPatches_[1] = std::make_shared<WeakeningRegion>(cuboidGeometries_[2]->weakeningRegions()[0]);
+ couplingIndices[1] = {2, 0};
+
+ nonmortarPatch_[2] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->right);
+ mortarPatch_[2] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->left);
+ weakPatches_[2] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[1]);
+ couplingIndices[2] = {1, 2};
+
+ for (size_t i=0; i<this->couplings_.size(); i++) {
+ auto& coupling = this->couplings_[i];
+ coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+ coupling->set(couplingIndices[i][0], couplingIndices[i][1], nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+ coupling->setWeakeningPatch(weakPatches_[i]);
+ coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale));
+ }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
+ using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+ std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
+
+ const double lengthScale = CuboidGeometry::lengthScale;
+
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& body = this->contactNetwork_.body(i);
+ const auto& leafVertexCount = body->nVertices();
+
+ // Neumann boundary
+ std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
+ body->addBoundaryCondition(neumannBoundary);
+
+ // upper Dirichlet Boundary
+ // identify Dirichlet nodes on leaf level
+ if (i>0) {
+ std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
+ for (int j=0; j<leafVertexCount; j++) {
+ if (leafFaces_[i]->upper.containsVertex(j))
+ (*velocityDirichletNodes)[j][0] = true;
+
+ #if MY_DIM == 3 //TODO: wrong, needs revision
+ if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
+ zeroDirichletNodes->at(j)[2] = true;
+ #endif
+ }
+
+ std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+
+ velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
+ velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+ body->addBoundaryCondition(velocityDirichletBoundary);
+ }
+ }
+
+ // lower Dirichlet Boundary
+ const auto& firstBody = this->contactNetwork_.body(0);
+ const auto& firstLeafVertexCount = firstBody->nVertices();
+ std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
+ for (int j=0; j<firstLeafVertexCount; j++) {
+ if (leafFaces_[0]->lower.containsVertex(j)) {
+ for (size_t d=0; d<dim; d++) {
+ (*zeroDirichletNodes)[j][d] = true;
+ }
+ }
+ }
+
+ std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
+
+ std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+ zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+ firstBody->addBoundaryCondition(zeroDirichletBoundary);
+}
+
+#include "threeblocksfactory_tmpl.cc"
diff --git a/src/factories/threeblocksfactory.hh b/src/factories/threeblocksfactory.hh
new file mode 100644
index 00000000..e3405572
--- /dev/null
+++ b/src/factories/threeblocksfactory.hh
@@ -0,0 +1,84 @@
+#ifndef SRC_THREEBLOCKSFACTORY_HH
+#define SRC_THREEBLOCKSFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../multi-body-problem-data/mybody.hh"
+#include "../multi-body-problem-data/grid/mygrids.hh"
+#include "../multi-body-problem-data/grid/cuboidgeometry.hh"
+
+template <class HostGridType, class VectorType>
+class ThreeBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+ using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+ using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+ using GlobalCoords = typename ContactNetwork::LocalVector;
+
+ using LeafGridView = typename ContactNetwork::GridView;
+ using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+ using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+ using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+ using LeafFaces = MyFaces<LeafGridView>;
+ using LevelFaces = MyFaces<LevelGridView>;
+
+ using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
+ using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
+
+ static const int dim = ContactNetwork::dim;
+
+ const std::shared_ptr<MyBodyData<dim>> bodyData_; // material properties of bodies
+
+ std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
+
+ std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
+
+ std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+ std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+ std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+ std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+ std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+ ThreeBlocksFactory(const Dune::ParameterTree& parset) :
+ Base(parset, 3, 3),
+ bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_, zenith_())),
+ cuboidGeometries_(3),
+ leafFaces_(3),
+ levelFaces_(3),
+ weakPatches_(3),
+ nonmortarPatch_(3),
+ mortarPatch_(3)
+ {}
+
+ void setBodies();
+ void setLevelBodies();
+ void setCouplings();
+ void setLevelCouplings() {}
+ void setBoundaryConditions();
+
+private:
+ static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+ #if MY_DIM == 2
+ return {0, 1};
+ #elif MY_DIM == 3
+ return {0, 1, 0};
+ #endif
+ }
+};
+#endif
+
+
diff --git a/src/factories/threeblocksfactory_tmpl.cc b/src/factories/threeblocksfactory_tmpl.cc
new file mode 100644
index 00000000..9ba11d9d
--- /dev/null
+++ b/src/factories/threeblocksfactory_tmpl.cc
@@ -0,0 +1,8 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class ThreeBlocksFactory<Grid, Vector>;
diff --git a/src/multi-body-problem-data/grid/cuboidgeometry.cc b/src/multi-body-problem-data/grid/cuboidgeometry.cc
index 77cd62f9..2fe7b2b5 100644
--- a/src/multi-body-problem-data/grid/cuboidgeometry.cc
+++ b/src/multi-body-problem-data/grid/cuboidgeometry.cc
@@ -12,76 +12,85 @@
#include "cuboidgeometry.hh"
-/*
-const CuboidGeometry::LocalVector CuboidGeometry::lift(const double v0, const double v1) {
- vc = 0;
- vc[0] = vc2D[0];
- vc[1] = vc2D[1];
-}
-
-const CuboidGeometry::LocalVector& CuboidGeometry::A() const { return A_; }
-const CuboidGeometry::LocalVector& CuboidGeometry::B() const { return B_; }
-const CuboidGeometry::LocalVector& CuboidGeometry::C() const { return C_; }
-const CuboidGeometry::LocalVector& CuboidGeometry::D() const { return D_; }
-const CuboidGeometry::LocalVector& CuboidGeometry::X() const { return X_; }
-const CuboidGeometry::LocalVector& CuboidGeometry::Y() const { return Y_; }
-double CuboidGeometry::depth() const { return depth_; }
-
-void CuboidGeometry::setupWeak(const LocalVector& weakOrigin, const double weakLength) {
- lift(weakOrigin, X_);
- const LocalVector2D Y({X_[0]+weakLength, X_[1]});
- lift(Y, Y_);
-}
-*/
-
#if MY_DIM == 3
-CuboidGeometry::CuboidGeometry(const LocalVector& origin,
- const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
- const double lowerWeakLength, const double upperWeakLength,
- const double length const double width, const double depth):
+template <class ctype>
+CuboidGeometry<ctype>::CuboidGeometry(const GlobalCoords& origin,
+ const double length = 1.00, const double height = 0.27, const double depth = 0.60) :
length_(length*lengthScale),
- width_(width*lengthScale),
- A(origin),
- B({origin[0]+length_, origin[1], 0}),
- C({origin[0]+length_, origin[1]+width_, 0}),
- D({origin[0], origin[1]+width_, 0}),
- V(lowerWeakOrigin),
- W({V[0]+lowerWeakLength*lengthScale, V[1], 0}),
- X(upperWeakOrigin),
- Y({X[0]+upperWeakLength*lengthScale, X[1], 0}),
- depth(depth_*lengthScale) {}
+ height_(height*lengthScale),
+ depth_(depth*lengthScale),
+ lowerLeft_(origin),
+ lowerRight_({origin[0]+length_, origin[1], 0}),
+ upperRight_({origin[0]+length_, origin[1]+height_, 0}),
+ upperLeft_({origin[0], origin[1]+height_, 0})
+ {}
#else
-CuboidGeometry::CuboidGeometry(const LocalVector& origin,
- const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
- const double lowerWeakLength, const double upperWeakLength,
- const double length, const double width):
+
+template <class ctype>
+CuboidGeometry<ctype>::CuboidGeometry(const GlobalCoords& origin,
+ const double length, const double height) :
length_(length*lengthScale),
- width_(width*lengthScale),
- A(origin),
- B({origin[0]+length_, origin[1]}),
- C({origin[0]+length_, origin[1]+width_}),
- D({origin[0], origin[1]+width_}),
- V(lowerWeakOrigin),
- W({V[0]+lowerWeakLength*lengthScale, V[1]}),
- X(upperWeakOrigin),
- Y({X[0]+upperWeakLength*lengthScale, X[1]}) {}
+ height_(height*lengthScale),
+ lowerLeft_(origin),
+ lowerRight_({origin[0]+length_, origin[1]}),
+ upperRight_({origin[0]+length_, origin[1]+height_}),
+ upperLeft_({origin[0], origin[1]+height_})
+ {}
#endif
+template <class ctype>
+void CuboidGeometry<ctype>::addWeakeningRegion(const WeakeningRegion& weakeningRegion) {
+ weakeningRegions_.emplace_back(weakeningRegion);
+}
+
+template <class ctype>
+void CuboidGeometry<ctype>::addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1) {
+ WeakeningRegion weakPatch;
+
+#if MY_DIM == 3 // TODO: Does not work yet
+ if (vertex0 != vertex1) {
+ weakPatch.vertices.resize(4);
+ weakPatch.vertices[0] = weakPatch.vertices[2] = vertex0;
+ weakPatch.vertices[1] = weakPatch.vertices[3] = vertex1;
+
+ for (size_t k = 0; k < 2; ++k) {
+ weakPatch.vertices[k][2] = -depth_ / 2.0;
+ weakPatch.vertices[k + 2][2] = depth_ / 2.0;
+ }
+ switch (parset.get<Config::PatchType>("patchType")) {
+ case Config::Rectangular:
+ break;
+ case Config::Trapezoidal:
+ weakPatch.vertices[1][0] += 0.05 * lengthScale;
+ weakPatch.vertices[3][0] -= 0.05 * lengthScale;
+ break;
+ default:
+ assert(false);
+ }
+ addWeakeningRegion(weakPatch);
+ }
+#else
+ if (vertex0 != vertex1) {
+ weakPatch.vertices.resize(2);
+ weakPatch.vertices[0] = vertex0;
+ weakPatch.vertices[1] = vertex1;
+ addWeakeningRegion(weakPatch);
+ }
+#endif
+}
-void CuboidGeometry::write() const {
- std::fstream writer("geometry", std::fstream::out);
- writer << "A = " << A << std::endl;
- writer << "B = " << B << std::endl;
- writer << "C = " << C << std::endl;
- writer << "D = " << D << std::endl;
- writer << "V = " << V << std::endl;
- writer << "W = " << W << std::endl;
- writer << "X = " << X << std::endl;
- writer << "Y = " << Y << std::endl;
+template <class ctype>
+void CuboidGeometry<ctype>::write() const {
+ std::fstream writer("geometry", std::fstream::out);
+ writer << "lowerLeft = " << lowerLeft_ << std::endl;
+ writer << "lowerRight = " << lowerRight_ << std::endl;
+ writer << "upperRight = " << upperRight_ << std::endl;
+ writer << "upperLeft = " << upperLeft_ << std::endl;
}
/*
-void CuboidGeometry::render() const {
+template <class ctype>
+void CuboidGeometry<ctype>::render() const {
#ifdef HAVE_CAIROMM
std::string const filename = "geometry.png";
double const width = 600;
@@ -202,3 +211,5 @@ void CuboidGeometry::render() const {
#endif
}
*/
+
+#include "cuboidgeometry_tmpl.cc"
diff --git a/src/multi-body-problem-data/grid/cuboidgeometry.hh b/src/multi-body-problem-data/grid/cuboidgeometry.hh
index 834f3f4d..56998355 100644
--- a/src/multi-body-problem-data/grid/cuboidgeometry.hh
+++ b/src/multi-body-problem-data/grid/cuboidgeometry.hh
@@ -1,51 +1,73 @@
#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
#define SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+#include <vector>
+
#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
-#include "../midpoint.hh"
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+template <class ctype = double>
class CuboidGeometry {
- public:
- typedef Dune::FieldVector<double, MY_DIM> LocalVector;
+public:
+ typedef Dune::FieldVector<ctype, MY_DIM> GlobalCoords;
+ using WeakeningRegion = ConvexPolyhedron<GlobalCoords>;
constexpr static double const lengthScale = 1.0; // scaling factor
- private:
- const double length_;
- const double width_;
- // const double lowerWeakLength_; // length of the lower weak zone
- // const double upperWeakLength_; // length of the lower weak zone
+private:
+ const ctype length_;
+ const ctype height_;
+#if MY_DIM == 3
+ const ctype depth_;
+#endif
- public:
// corners of cube
- const LocalVector A;
- const LocalVector B;
- const LocalVector C;
- const LocalVector D;
+ const GlobalCoords lowerLeft_;
+ const GlobalCoords lowerRight_;
+ const GlobalCoords upperRight_;
+ const GlobalCoords upperLeft_;
- // corners of lower weakening patch
- const LocalVector V;
- const LocalVector W;
-
- // corners of upper weakening patch
- const LocalVector X;
- const LocalVector Y;
+ // weakening regions
+ std::vector<WeakeningRegion> weakeningRegions_;
+public:
#if MY_DIM == 3
- const double depth;
+ CuboidGeometry(const GlobalCoords& origin,
+ const double length = 1.00, const double height = 0.27, const double depth = 0.60);
- CuboidGeometry(const LocalVector& origin,
- const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
- const double lowerWeakLength = 0.20, const double upperWeakLength = 0.20,
- const double length = 1.00, const double width = 0.27, const double depth = 0.60);
+ const auto& depth() const {
+ return depth_;
+ }
#else
- CuboidGeometry(const LocalVector& origin,
- const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
- const double lowerWeakLength = 0.20, const double upperWeakLength = 0.20,
- const double length = 1.00, const double width = 0.27);
+ CuboidGeometry(const GlobalCoords& origin,
+ const double length = 1.00, const double height = 0.27);
#endif
+ void addWeakeningRegion(const WeakeningRegion& weakeningRegion);
+ void addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1);
+
+ const auto& lowerLeft() const {
+ return lowerLeft_;
+ }
+
+ const auto& lowerRight() const {
+ return lowerRight_;
+ }
+
+ const auto& upperRight() const {
+ return upperRight_;
+ }
+
+ const auto& upperLeft() const {
+ return upperLeft_;
+ }
+
+ const auto& weakeningRegions() const {
+ return weakeningRegions_;
+ }
+
void write() const;
// void render() const;
diff --git a/src/multi-body-problem-data/grid/cuboidgeometry_tmpl.cc b/src/multi-body-problem-data/grid/cuboidgeometry_tmpl.cc
new file mode 100644
index 00000000..6065a02b
--- /dev/null
+++ b/src/multi-body-problem-data/grid/cuboidgeometry_tmpl.cc
@@ -0,0 +1,8 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "cuboidgeometry.hh"
+
+template class CuboidGeometry<typename Grid::ctype>;
diff --git a/src/multi-body-problem-data/grid/mygrids.cc b/src/multi-body-problem-data/grid/mygrids.cc
index 3d2fafda..e18f316b 100644
--- a/src/multi-body-problem-data/grid/mygrids.cc
+++ b/src/multi-body-problem-data/grid/mygrids.cc
@@ -11,7 +11,7 @@
#include "simplexmanager.hh"
// Fix: 3D case (still Elias code)
-template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_) :
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_) :
cuboidGeometries(cuboidGeometries_)
{
size_t const gridCount = cuboidGeometries.size();
@@ -21,10 +21,10 @@ template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::
for (size_t idx=0; idx<grids.size(); idx++) {
const auto& cuboidGeometry = *cuboidGeometries[idx];
- const auto& A = cuboidGeometry.A;
- const auto& B = cuboidGeometry.B;
- const auto& C = cuboidGeometry.C;
- const auto& D = cuboidGeometry.D;
+ const auto& A = cuboidGeometry.lowerLeft();
+ const auto& B = cuboidGeometry.lowerRight();
+ const auto& C = cuboidGeometry.upperRight();
+ const auto& D = cuboidGeometry.upperLeft();
unsigned int const vc = 4;
@@ -51,8 +51,8 @@ template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::
#if MY_DIM == 3
for (size_t k = 0; k < vc; ++k) {
- vertices[k][2] = -cuboidGeometry.depth / 2.0;
- vertices[k + vc][2] = cuboidGeometry.depth / 2.0;
+ vertices[k][2] = -cuboidGeometry.depth() / 2.0;
+ vertices[k + vc][2] = cuboidGeometry.depth() / 2.0;
}
#endif
@@ -96,7 +96,7 @@ std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() {
template <class Grid>
template <class GridView>
MyFaces<GridView> GridsConstructor<Grid>::constructFaces(
- GridView const &gridView, CuboidGeometry const &cuboidGeometry) {
+ GridView const &gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry) {
return MyFaces<GridView>(gridView, cuboidGeometry);
}
@@ -132,7 +132,7 @@ bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
}
template <class GridView>
-MyFaces<GridView>::MyFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry_)
+MyFaces<GridView>::MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_)
:
#if MY_DIM == 3
lower(gridView),
@@ -150,28 +150,28 @@ MyFaces<GridView>::MyFaces(GridView const &gridView, CuboidGeometry const &cuboi
cuboidGeometry(cuboidGeometry_)
{
lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return xyBetween(cuboidGeometry.A, cuboidGeometry.B, in.geometry().center());
+ return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.lowerRight(), in.geometry().center());
});
right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return xyBetween(cuboidGeometry.B, cuboidGeometry.C, in.geometry().center());
+ return xyBetween(cuboidGeometry.lowerRight(), cuboidGeometry.upperRight(), in.geometry().center());
});
upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return xyBetween(cuboidGeometry.D, cuboidGeometry.C, in.geometry().center());
+ return xyBetween(cuboidGeometry.upperLeft(), cuboidGeometry.upperRight(), in.geometry().center());
});
left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return xyBetween(cuboidGeometry.A, cuboidGeometry.D, in.geometry().center());
+ return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.upperLeft(), in.geometry().center());
});
#if MY_DIM == 3
front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return isClose(cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+ 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]);
+ return isClose(-cuboidGeometry.depth() / 2.0, in.geometry().center()[2]);
});
#endif
}
diff --git a/src/multi-body-problem-data/grid/mygrids.hh b/src/multi-body-problem-data/grid/mygrids.hh
index f37716eb..574efa1d 100644
--- a/src/multi-body-problem-data/grid/mygrids.hh
+++ b/src/multi-body-problem-data/grid/mygrids.hh
@@ -20,10 +20,10 @@ template <class GridView> struct MyFaces {
BoundaryPatch<GridView> back;
#endif
- MyFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry_);
+ MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_);
private:
- CuboidGeometry const &cuboidGeometry;
+ const CuboidGeometry<typename GridView::ctype>& cuboidGeometry;
bool isClose(double a, double b) {
return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale;
@@ -46,15 +46,15 @@ template <class GridView> struct MyFaces {
template <class Grid> class GridsConstructor {
public:
- GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
+ GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_);
std::vector<std::shared_ptr<Grid>>& getGrids();
template <class GridView>
- MyFaces<GridView> constructFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry);
+ MyFaces<GridView> constructFaces(const GridView& gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry);
private:
- std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries;
+ const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries;
std::vector<Dune::GridFactory<Grid>> gridFactories;
std::vector<std::shared_ptr<Grid>> grids;
};
diff --git a/src/multi-body-problem-data/grid/mygrids_tmpl.cc b/src/multi-body-problem-data/grid/mygrids_tmpl.cc
index 32a19759..e626dd27 100644
--- a/src/multi-body-problem-data/grid/mygrids_tmpl.cc
+++ b/src/multi-body-problem-data/grid/mygrids_tmpl.cc
@@ -12,10 +12,10 @@ template struct MyFaces<DeformedGrid::LeafGridView>;
template struct MyFaces<DeformedGrid::LevelGridView>;
template MyFaces<DeformedGrid::LeafGridView> GridsConstructor<Grid>::constructFaces(
- DeformedGrid::LeafGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+ DeformedGrid::LeafGridView const &gridView, CuboidGeometry<typename Grid::ctype> const &CuboidGeometry_);
template MyFaces<DeformedGrid::LevelGridView> GridsConstructor<Grid>::constructFaces(
- DeformedGrid::LevelGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+ DeformedGrid::LevelGridView const &gridView, CuboidGeometry<typename Grid::ctype> const &CuboidGeometry_);
template void refine<Grid, LocalVector>(
Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
diff --git a/src/multi-body-problem-data/weakpatch.hh b/src/multi-body-problem-data/weakpatch.hh
deleted file mode 100644
index 047dc482..00000000
--- a/src/multi-body-problem-data/weakpatch.hh
+++ /dev/null
@@ -1,67 +0,0 @@
-#ifndef SRC_MULTI_BODY_PROBLEM_DATA_WEAKPATCH_HH
-#define SRC_MULTI_BODY_PROBLEM_DATA_WEAKPATCH_HH
-
-#include <dune/common/parametertree.hh>
-
-#include <dune/fufem/geometry/convexpolyhedron.hh>
-
-#include "grid/cuboidgeometry.hh"
-
-template <class LocalVector>
-void getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry, ConvexPolyhedron<LocalVector>& lowerWeakPatch, ConvexPolyhedron<LocalVector>& upperWeakPatch) {
-
-#if MY_DIM == 3 // TODO: Does not work yet
- if (cuboidGeometry.V != cuboidGeometry.W) {
- lowerWeakPatch.vertices.resize(4);
- lowerWeakPatch.vertices[0] = lowerWeakPatch.vertices[2] = cuboidGeometry.V;
- lowerWeakPatch.vertices[1] = lowerWeakPatch.vertices[3] = cuboidGeometry.W;
- for (size_t k = 0; k < 2; ++k) {
- lowerWeakPatch.vertices[k][2] = -cuboidGeometry.depth / 2.0;
- lowerWeakPatch.vertices[k + 2][2] = cuboidGeometry.depth / 2.0;
- }
- switch (parset.get<Config::PatchType>("patchType")) {
- case Config::Rectangular:
- break;
- case Config::Trapezoidal:
- lowerWeakPatch.vertices[1][0] += 0.05 * cuboidGeometry.lengthScale;
- lowerWeakPatch.vertices[3][0] -= 0.05 * cuboidGeometry.lengthScale;
- break;
- default:
- assert(false);
- }
- }
-
- if (cuboidGeometry.X != cuboidGeometry.Y) {
- upperWeakPatch.vertices.resize(4);
- upperWeakPatch.vertices[0] = upperWeakPatch.vertices[2] = cuboidGeometry.X;
- upperWeakPatch.vertices[1] = upperWeakPatch.vertices[3] = cuboidGeometry.Y;
- for (size_t k = 0; k < 2; ++k) {
- upperWeakPatch.vertices[k][2] = -cuboidGeometry.depth / 2.0;
- upperWeakPatch.vertices[k + 2][2] = cuboidGeometry.depth / 2.0;
- }
- switch (parset.get<Config::PatchType>("patchType")) {
- case Config::Rectangular:
- break;
- case Config::Trapezoidal:
- upperWeakPatch.vertices[1][0] += 0.05 * cuboidGeometry.lengthScale;
- upperWeakPatch.vertices[3][0] -= 0.05 * cuboidGeometry.lengthScale;
- break;
- default:
- assert(false);
- }
- }
-#else
- if (cuboidGeometry.V != cuboidGeometry.W) {
- lowerWeakPatch.vertices.resize(2);
- lowerWeakPatch.vertices[0] = cuboidGeometry.V;
- lowerWeakPatch.vertices[1] = cuboidGeometry.W;
- }
-
- if (cuboidGeometry.X != cuboidGeometry.Y) {
- upperWeakPatch.vertices.resize(2);
- upperWeakPatch.vertices[0] = cuboidGeometry.X;
- upperWeakPatch.vertices[1] = cuboidGeometry.Y;
- }
-#endif
-};
-#endif
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index e20cc637..6ab6c3dd 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -55,11 +55,11 @@
#include "data-structures/enumparser.hh"
#include "data-structures/enums.hh"
-#include "data-structures/levelcontactnetwork.hh"
+#include "data-structures/contactnetwork.hh"
#include "data-structures/matrices.hh"
#include "data-structures/program_state.hh"
-#include "factories/stackedblocksfactory.hh"
+#include "factories/threeblocksfactory.hh"
//#include "io/hdf5-levelwriter.hh"
#include "io/hdf5/restart-io.hh"
@@ -70,6 +70,7 @@
#include "multi-body-problem-data/grid/mygrids.hh"
#include "spatial-solving/tnnmg/functional.hh"
+#include "spatial-solving/preconditioners/supportpatchfactory.hh"
#include "spatial-solving/tnnmg/localbisectionsolver.hh"
#include "spatial-solving/contact/nbodyassembler.hh"
#include "spatial-solving/solverfactory.hh"
@@ -90,6 +91,34 @@
size_t const dims = MY_DIM;
+template <class SupportPatchFactory>
+void testSuite(const SupportPatchFactory& supportPatchFactory, const size_t bodyID, const int patchDepth = 0) {
+/* const auto& coarseContactNetwork = supportPatchFactory.coarseContactNetwork();
+ const auto& gridView = coarseContactNetwork.body(bodyID)->gridView();
+
+ for (const auto& e : elements(gridView)) {
+
+ }
+
+ using Patch = typename SupportPatchFactory::Patch;
+ Patch patch0, patch1, patch2, patch3;
+
+ // (const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0)
+
+ // interior patch inside of one body
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch0, patchDepth);
+
+ // patch at friction interface with two bodies in contact
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch1 patchDepth);
+
+ // patch at friction interface with two bodies in contact and with dirichlet boundary
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch2, patchDepth);
+
+ // crosspoint patch, where all 3 bodies are in contact
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch3, patchDepth);*/
+}
+
+
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
@@ -126,41 +155,99 @@ int main(int argc, char *argv[]) {
// ----------------------
// set up contact network
// ----------------------
- StackedBlocksFactory<Grid, dims> stackedBlocksFactory(parset);
- using LevelContactNetwork = typename StackedBlocksFactory<Grid, dims>::LevelContactNetwork;
- stackedBlocksFactory.build();
+ ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
+ using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
+ threeBlocksFactory.build();
- LevelContactNetwork& levelContactNetwork = stackedBlocksFactory.levelContactNetwork();
- const size_t bodyCount = levelContactNetwork.nBodies();
+ ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork();
+ const size_t bodyCount = contactNetwork.nBodies();
for (size_t i=0; i<bodyCount; i++) {
- printDofLocation(levelContactNetwork.leafView(i));
- /* Vector def(levelContactNetwork.deformedGrids()[i]->size(dims));
+ // printDofLocation(contactNetwork.body(i)->gridView());
+ /* Vector def(contactNetwork.deformedGrids()[i]->size(dims));
def = 1;
deformedGridComplex.setDeformation(def, i);*/
- /*auto& levelViews = levelContactNetwork.levelViews(i);
+ /*auto& levelViews = contactNetwork.levelViews(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"); */
+ writeToVTK(contactNetwork.leafView(i), "", "body_" + std::to_string(i) + "_leaf"); */
}
// ----------------------------
- // assemble levelContactNetwork
+ // assemble contactNetwork
// ----------------------------
+ contactNetwork.assemble();
+
+ const auto & coarseContactNetwork = *contactNetwork.level(0);
+ const auto & fineContactNetwork = *contactNetwork.level(1);
+ SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
+ size_t bodyID = 0;
+ size_t patchDepth = 0;
+
+ writeToVTK(fineContactNetwork.body(bodyID)->gridView(), "", "body_" + std::to_string(bodyID) + "_fine");
+ writeToVTK(coarseContactNetwork.body(bodyID)->gridView(), "", "body_" + std::to_string(bodyID) + "_coarse");
+
+ using Patch = typename SupportPatchFactory<typename ContactNetwork::LevelContactNetwork>::Patch;
+ Patch patch;
+
+ const auto& gridView = coarseContactNetwork.body(bodyID)->gridView();
+
+ Dune::PQkLocalFiniteElementCache<double, double, dims, 1> cache;
+ Dune::BitSetVector<1> vertexVisited(gridView.size(dims));
+ vertexVisited.unsetAll();
- levelContactNetwork.assemble();
- //printMortarBasis<Vector>(levelContactNetwork.nBodyAssembler());
+ for (const auto& e: elements(gridView)) {
+ const auto& refElement = Dune::ReferenceElements<double, dims>::general(e.type());
+
+ for (size_t i=0; i<refElement.size(dims); i++) {
+ auto localIdx = cache.get(e.type()).localCoefficients().localKey(i).subEntity();
+ auto globalIdx = gridView.indexSet().subIndex(e, i, dims);
+
+ if (!vertexVisited[globalIdx][0]) {
+ vertexVisited[globalIdx][0] = true;
+ supportPatchFactory.build(bodyID, e, i, patch, patchDepth);
+
+ print(patch, "patch:");
+
+ size_t c = 0;
+ for (size_t j=0; j<bodyCount; j++) {
+ const auto& gv = fineContactNetwork.body(j)->gridView();
+
+ ScalarVector patchVec(gv.size(dims));
+ for (size_t l=0; l<patchVec.size(); l++) {
+ if (patch[c++][0]) {
+ patchVec[l][0] = 1;
+ }
+ }
+
+ print(patchVec, "patchVec");
+
+ // output patch
+ writeToVTK(gv, patchVec, "", "patch_" + std::to_string(globalIdx) + "_body_" + std::to_string(j));
+ }
+ }
+
+
+ }
+ }
+ return 1;
+
+ //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
// -----------------
// init input/output
// -----------------
- const auto& leafVertexCounts = levelContactNetwork.leafVertexCounts();
+ std::vector<size_t> nVertices(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ nVertices[i] = contactNetwork.body(i)->nVertices();
+ }
+
using MyProgramState = ProgramState<Vector, ScalarVector>;
- MyProgramState programState(leafVertexCounts);
+ MyProgramState programState(nVertices);
auto const firstRestart = parset.get<size_t>("io.restarts.first");
auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
auto const writeRestarts = parset.get<bool>("io.restarts.write");
@@ -181,18 +268,18 @@ int main(int argc, char *argv[]) {
auto restartIO = handleRestarts
? std::make_unique<RestartIO<MyProgramState>>(
- *restartFile, leafVertexCounts)
+ *restartFile, nVertices)
: nullptr;
if (firstRestart > 0) // automatically adjusts the time and timestep
restartIO->read(firstRestart, programState);
else
- programState.setupInitialConditions(parset, levelContactNetwork);
+ programState.setupInitialConditions(parset, contactNetwork);
- auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
+ auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
- levelContactNetwork.body(i)->setDeformation(programState.u[i]);
+ contactNetwork.body(i)->setDeformation(programState.u[i]);
}
nBodyAssembler.assembleTransferOperator();
nBodyAssembler.assembleObstacle();
@@ -200,9 +287,9 @@ int main(int argc, char *argv[]) {
// ------------------------
// assemble global friction
// ------------------------
- levelContactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+ contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
- auto& globalFriction = levelContactNetwork.globalFriction();
+ auto& globalFriction = contactNetwork.globalFriction();
globalFriction.updateAlpha(programState.alpha);
using MyVertexBasis = typename Assembler::VertexBasis;
@@ -212,21 +299,21 @@ int main(int argc, char *argv[]) {
std::vector<const MyCellBasis* > cellBases(bodyCount);
for (size_t i=0; i<bodyCount; i++) {
- const auto& body = levelContactNetwork.body(i);
+ const auto& body = contactNetwork.body(i);
vertexBases[i] = &(body->assembler()->vertexBasis);
cellBases[i] = &(body->assembler()->cellBasis);
auto& vertexCoords = vertexCoordinates[i];
- vertexCoords.resize(leafVertexCounts[i]);
+ vertexCoords.resize(nVertices[i]);
Dune::MultipleCodimMultipleGeomTypeMapper<
- DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->leafView(), Dune::mcmgVertexLayout());
- for (auto &&v : vertices(body->leafView()))
+ DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(body->gridView()))
vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
- //typename LevelContactNetwork::BoundaryPatches frictionBoundaries;
- //levelContactNetwork.boundaryPatches("friction", frictionBoundaries);
+ //typename contactNetwork::BoundaryPatches frictionBoundaries;
+ //contactNetwork.boundaryPatches("friction", frictionBoundaries);
/*
auto dataWriter =
@@ -242,7 +329,7 @@ int main(int argc, char *argv[]) {
auto const report = [&](bool initial = false) {
/*if (writeData) {
- dataWriter->reportSolution(programState, levelContactNetwork.globalFriction());
+ dataWriter->reportSolution(programState, contactNetwork.globalFriction());
if (!initial)
dataWriter->reportIterations(programState, iterationCount);
dataFile->flush();
@@ -264,7 +351,7 @@ int main(int argc, char *argv[]) {
std::vector<ScalarVector> stress(bodyCount);
for (size_t i=0; i<bodyCount; i++) {
- const auto& body = levelContactNetwork.body(i);
+ const auto& body = contactNetwork.body(i);
body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
body->data()->getPoissonRatio(),
programState.u[i], stress[i]);
@@ -282,7 +369,7 @@ int main(int argc, char *argv[]) {
// -------------------
BitVector totalDirichletNodes;
- levelContactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+ contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
/*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
bool val = false;
@@ -301,16 +388,16 @@ int main(int argc, char *argv[]) {
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
using NonlinearFactory = SolverFactory<Functional, BitVector>;
- using BoundaryFunctions = typename LevelContactNetwork::BoundaryFunctions;
- using BoundaryNodes = typename LevelContactNetwork::BoundaryNodes;
+ using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+ using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
StateUpdater<ScalarVector, Vector>>;
BoundaryFunctions velocityDirichletFunctions;
- levelContactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+ contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
BoundaryNodes dirichletNodes;
- levelContactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+ contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
@@ -319,7 +406,7 @@ int main(int argc, char *argv[]) {
}
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
- levelContactNetwork.frictionNodes(frictionNodes);
+ contactNetwork.frictionNodes(frictionNodes);
for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
@@ -330,7 +417,7 @@ int main(int argc, char *argv[]) {
parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunctions,
dirichletNodes,
- levelContactNetwork.matrices(),
+ contactNetwork.matrices(),
programState.u,
programState.v,
programState.a),
@@ -338,13 +425,13 @@ int main(int argc, char *argv[]) {
parset.get<Config::stateModel>("boundary.friction.stateModel"),
programState.alpha,
nBodyAssembler.getContactCouplings(),
- levelContactNetwork.couplings())
+ contactNetwork.couplings())
);
auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
- const auto& stateEnergyNorms = levelContactNetwork.stateEnergyNorms();
+ const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
auto const mustRefine = [&](Updaters &coarseUpdater,
Updaters &fineUpdater) {
@@ -387,10 +474,10 @@ int main(int argc, char *argv[]) {
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
- typename LevelContactNetwork::ExternalForces externalForces;
- levelContactNetwork.externalForces(externalForces);
+ typename ContactNetwork::ExternalForces externalForces;
+ contactNetwork.externalForces(externalForces);
- AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(nBodyAssembler)>, Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
+ AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(nBodyAssembler)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
adaptiveTimeStepper(parset, nBodyAssembler, totalDirichletNodes, globalFriction, frictionNodes, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
@@ -412,11 +499,7 @@ int main(int argc, char *argv[]) {
print(programState.a, "current a:");
print(programState.alpha, "current alpha:");
- for (size_t i=0; i<bodyCount; i++) {
- levelContactNetwork.body(i)->setDeformation(programState.u[i]);
- }
- nBodyAssembler.assembleTransferOperator();
- nBodyAssembler.assembleObstacle();
+ contactNetwork.setDeformation(programState.u);
report();
diff --git a/src/spatial-solving/contact/dualmortarcoupling.hh b/src/spatial-solving/contact/dualmortarcoupling.hh
index 90bc49dc..81c6e628 100644
--- a/src/spatial-solving/contact/dualmortarcoupling.hh
+++ b/src/spatial-solving/contact/dualmortarcoupling.hh
@@ -93,12 +93,12 @@ class DualMortarCoupling {
}
/** \brief Setup leaf nonmortar and mortar patches. */
- /*virtual void setupContactPatch(const LevelBoundaryPatch0& coarseNonmortar,
+ virtual void setupContactPatch(const LevelBoundaryPatch0& coarseNonmortar,
const LevelBoundaryPatch1& coarseMortar) {
BoundaryPatchProlongator<GridType0>::prolong(coarseNonmortar,nonmortarBoundary_);
BoundaryPatchProlongator<GridType1>::prolong(coarseMortar,mortarBoundary_);
- }*/
+ }
/** \brief Return the nonmortar boundary. */
const BoundaryPatch0& nonmortarBoundary() const {
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index f600e878..2c741b77 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -37,14 +37,14 @@ void FixedPointIterationCounter::operator+=(
multigridIterations += other.multigridIterations;
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIterator(
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIterator(
Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
- const std::vector<const ErrorNorm* >& errorNorms)
+ const ErrorNorms& errorNorms)
: parset_(parset),
nBodyAssembler_(nBodyAssembler),
ignoreNodes_(ignoreNodes),
@@ -58,9 +58,9 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIter
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
errorNorms_(errorNorms) {}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
FixedPointIterationCounter
-FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
+FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
@@ -237,8 +237,8 @@ std::ostream &operator<<(std::ostream &stream,
<< ")";
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeVelocities(
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::relativeVelocities(
const Vector& v,
std::vector<Vector>& v_rel) const {
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index fa93c087..867e8058 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -23,13 +23,13 @@ struct FixedPointIterationCounter {
std::ostream &operator<<(std::ostream &stream,
FixedPointIterationCounter const &fpic);
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
class FixedPointIterator {
using Functional = typename Factory::Functional;
using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
- using Nonlinearity = typename Factory::Functional::Nonlinearity;
+ using Nonlinearity = typename Factory::Functional::Nonlinearity;
const static int dims = Vector::block_type::dimension;
@@ -50,7 +50,7 @@ class FixedPointIterator {
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
- const std::vector<const ErrorNorm* >& errorNorms);
+ const ErrorNorms& errorNorms);
FixedPointIterationCounter run(Updaters updaters,
const std::vector<Matrix>& velocityMatrices,
@@ -72,6 +72,6 @@ class FixedPointIterator {
size_t velocityMaxIterations_;
double velocityTolerance_;
Solver::VerbosityMode verbosity_;
- const std::vector<const ErrorNorm* >& errorNorms_;
+ const ErrorNorms& errorNorms_;
};
#endif
diff --git a/src/spatial-solving/fixedpointiterator_tmpl.cc b/src/spatial-solving/fixedpointiterator_tmpl.cc
index a91a7c6a..33acfce6 100644
--- a/src/spatial-solving/fixedpointiterator_tmpl.cc
+++ b/src/spatial-solving/fixedpointiterator_tmpl.cc
@@ -8,21 +8,21 @@
#include <dune/solvers/norms/energynorm.hh>
#include "../spatial-solving/solverfactory_tmpl.cc"
-#include "../data-structures/levelcontactnetwork_tmpl.cc"
+#include "../data-structures/contactnetwork_tmpl.cc"
#include "../time-stepping/rate/rateupdater.hh"
#include "../time-stepping/state/stateupdater.hh"
#include "../time-stepping/updaters.hh"
-using NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
-using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
-using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
+using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
+template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index 76d64e37..75db9ad2 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -7,6 +7,7 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/fvector.hh>
+#include <dune/common/hybridutilities.hh>
#include <dune/fufem/referenceelementhelper.hh>
#include <dune/grid/common/mcmgmapper.hh>
@@ -19,12 +20,13 @@
template <class LevelContactNetwork>
class NetworkIndexMapper {
private:
+ static const int dim = LevelContactNetwork::dim; //TODO
+ using ctype = typename LevelContactNetwork::GridView::ctype;
+
using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
- using GridType = typename LevelContactNetwork::DeformedGridType;
- using Element = typename GridType::template Codim<0>::Entity;
+ using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
- static const int dim = LevelContactNetwork::dimension; //TODO
- static const FiniteElementCache cache_;
+ const FiniteElementCache cache_;
const LevelContactNetwork& levelContactNetwork_;
const size_t nBodies_;
@@ -62,7 +64,7 @@ class NetworkIndexMapper {
size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
- auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
+ auto localIdx = cache().get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
}
@@ -82,11 +84,13 @@ class NetworkIndexMapper {
return elementOffSet_[bodyID] + gridView.indexSet().index(elem);
}
- static const auto& cache() const {
+ const auto& cache() const {
return cache_;
}
};
+
+
/**
* constructs BitSetVector whose entries are
* true, if vertex is outside of patch or part of patch boundary
@@ -99,6 +103,11 @@ class SupportPatchFactory
using Patch = Dune::BitSetVector<1>;
private:
+ static const int dim = LevelContactNetwork::dim; //TODO
+ using ctype = typename LevelContactNetwork::ctype;
+
+ using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
+
struct BodyElement {
size_t bodyID;
Element element;
@@ -111,40 +120,33 @@ class SupportPatchFactory
struct RemoteIntersection {
size_t bodyID; // used to store bodyID of outside elem
- GridGlue& glue;
+ size_t contactCouplingID;
size_t intersectionID;
bool flip = false;
- void set(const size_t _bodyID, const GridGlue& _glue, const size_t _intersectionID, bool _flip = false) {
+ void set(const size_t _bodyID, const size_t _contactCouplingID, const size_t _intersectionID, const bool _flip = false) {
bodyID = _bodyID;
- glue = _glue;
+ contactCouplingID = _contactCouplingID;
intersectionID = _intersectionID;
flip = _flip;
}
- auto& get() {
- return (flip) ? glue.getIntersection(intersectionID).flip() : glue.getIntersection(intersectionID);
+ template <class GridGlue>
+ const auto& get(const GridGlue& glue) const {
+ return glue.getIntersection(intersectionID);
}
};
- static const int dim = LevelContactNetwork::dimension; //TODO
- using ctype = typename LevelContactNetwork::ctype;
-
- //using Vertex = typename GridType::template Codim<dim>::Entity;
-
- using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity;
- using BodyElement = BodyElement<Element>;
-
using NetworkIndexMapper = NetworkIndexMapper<LevelContactNetwork>;
- using GridGlue = ;
-
const size_t nBodies_;
const LevelContactNetwork& coarseContactNetwork_;
const LevelContactNetwork& fineContactNetwork_;
+ size_t nVertices_ = 0;
+
std::map<size_t, std::vector<RemoteIntersection>> coarseIntersectionMapper_; // map faceID to neighbor elements
std::map<size_t, std::vector<RemoteIntersection>> fineIntersectionMapper_; // map faceID to neighbor elements
std::map<size_t, std::vector<size_t>> neighborFaceDofs_; // map faceID to neighbor element dofs
@@ -156,7 +158,7 @@ class SupportPatchFactory
public:
SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
- nBodies_(coarseContactNetwork_.nBodies()),
+ nBodies_(coarseContactNetwork.nBodies()),
coarseContactNetwork_(coarseContactNetwork),
fineContactNetwork_(fineContactNetwork),
coarseIndices_(coarseContactNetwork_),
@@ -164,14 +166,17 @@ class SupportPatchFactory
assert(nBodies_ == fineContactNetwork_.nBodies());
+ for (size_t i=0; i<nBodies_; i++) {
+ nVertices_ += fineContactNetwork_.body(i)->nVertices();
+ }
+
setup(coarseContactNetwork_, coarseIndices_, coarseIntersectionMapper_);
setup(fineContactNetwork_, fineIndices_, fineIntersectionMapper_);
// setup neighborFaceDofs_
- const auto& contactCouplings = coarseContactNetwork_.nBodyAssembler().getContactCouplings();
- for (size_t i=0; i<contactCouplings.size(); i++) {
- const auto& coupling = nBodyAssembler.getCoupling(i);
- const auto& glue = *contactCouplings[i].getGlue();
+ for (size_t i=0; i<coarseContactNetwork_.nCouplings(); i++) {
+ const auto& coupling = *coarseContactNetwork_.coupling(i);
+ const auto& glue = *coarseContactNetwork_.glue(i);
const size_t nmBody = coupling.gridIdx_[0];
const size_t mBody = coupling.gridIdx_[1];
@@ -180,7 +185,7 @@ class SupportPatchFactory
size_t nmFaceIdx = coarseIndices_.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
size_t mFaceIdx = coarseIndices_.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
- std::array<std::set<size_t>, 2> dofs;
+ std::vector<std::set<size_t>> dofs;
remoteIntersectionDofs(coarseIndices_, rIs, nmBody, mBody, dofs);
neighborFaceDofs_[nmFaceIdx].insert(neighborFaceDofs_[nmFaceIdx].end(), dofs[1].begin(), dofs[1].end());
neighborFaceDofs_[mFaceIdx].insert(neighborFaceDofs_[mFaceIdx].end(), dofs[0].begin(), dofs[0].end());
@@ -207,10 +212,16 @@ class SupportPatchFactory
*/
}
- void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0) {
- BodyElement seedElement(bodyID, coarseElement);
+ void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const size_t patchDepth = 0) {
+ std::cout << "Building SupportPatch... ";
+
+ patchDofs.resize(nVertices_);
+ patchDofs.setAll();
+
+ BodyElement seedElement;
+ seedElement.set(bodyID, coarseElement);
- auto isPatchIntersection = [](auto& is, const size_t bodyID, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
+ auto isPatchIntersection = [this](auto& is, const size_t bodyID, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
std::set<size_t> dofs;
intersectionDofs(coarseIndices_, is, bodyID, dofs);
@@ -255,10 +266,11 @@ class SupportPatchFactory
const auto& elem = patchSeed.element;
const auto& gridView = coarseContactNetwork_.body(bodyIdx)->gridView();
- for (auto& it : intersections(gridView, elem)) {
+ for (const auto& it : intersections(gridView, elem)) {
if (isPatchIntersection(it, bodyIdx, patchVertices, newPatchVertices)) {
if (it.neighbor()) {
- BodyElement neighbor(bodyIdx, it.outside());
+ BodyElement neighbor;
+ neighbor.set(bodyIdx, it.outside());
if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
@@ -272,12 +284,13 @@ class SupportPatchFactory
const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
BodyElement neighbor;
- neighbor.set(intersections[i].bodyID, intersections[i].get().outside());
+ const auto& glue = *coarseContactNetwork_.glue(intersections[i].contactCouplingID);
+ neighbor.set(intersections[i].bodyID, intersections[i].get(glue).outside());
if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
- patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
+ patchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
patchElements.emplace_back(neighbor);
patchSeeds.push(neighbor);
}
@@ -285,7 +298,8 @@ class SupportPatchFactory
}
} else {
if (it.neighbor()) {
- BodyElement neighbor(bodyIdx, it.outside());
+ BodyElement neighbor;
+ neighbor.set(bodyIdx, it.outside());
if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
@@ -298,12 +312,13 @@ class SupportPatchFactory
const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
BodyElement neighbor;
- neighbor.set(intersections[i].bodyID, intersections[i].get().outside());
+ const auto& glue = *coarseContactNetwork_.glue(intersections[i].contactCouplingID);
+ neighbor.set(intersections[i].bodyID, intersections[i].get(glue).outside());
if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
- newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
+ newPatchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
nextSeeds.emplace_back(neighbor);
}
}
@@ -319,13 +334,13 @@ class SupportPatchFactory
}
// construct fine patch
- using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::GridType>;
std::set<size_t> patchBoundary;
for (size_t i=0; i<patchElements.size(); ++i) {
const auto& coarseElem = patchElements[i];
const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
- const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->gridView().level();
+ const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->level();
// add fine patch elements
HierarchicLevelIteratorType endHierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
@@ -335,16 +350,19 @@ class SupportPatchFactory
addLocalDofs(coarseElem, fineElem, visitedElements, patchBoundary, patchDofs);
}
}
+
+ std::cout << "Success!" << std::endl;
+ }
+
+ const auto& coarseContactNetwork() {
+ return coarseContactNetwork_;
}
private:
void setup(const LevelContactNetwork& levelContactNetwork, const NetworkIndexMapper& indices, std::map<size_t, std::vector<RemoteIntersection>>& faceToRIntersections) {
- const auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
- const auto& contactCouplings = nBodyAssembler.getContactCouplings();
-
- for (size_t i=0; i<contactCouplings.size(); i++) {
- const auto& coupling = nBodyAssembler.getCoupling(i);
- const auto& glue = *contactCouplings[i].getGlue();
+ for (size_t i=0; i<levelContactNetwork.nCouplings(); i++) {
+ const auto& coupling = *levelContactNetwork.coupling(i);
+ const auto& glue = *levelContactNetwork.glue(i);
const size_t nmBody = coupling.gridIdx_[0];
const size_t mBody = coupling.gridIdx_[1];
@@ -354,10 +372,12 @@ class SupportPatchFactory
size_t nmFaceIdx = indices.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
size_t mFaceIdx = indices.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
- RemoteIntersection nmIntersection(mBody, glue, rIsID, false);
+ RemoteIntersection nmIntersection;
+ nmIntersection.set(mBody, i, rIsID, false);
faceToRIntersections[nmFaceIdx].emplace_back(nmIntersection);
- RemoteIntersection mIntersection(nmBody, glue, rIsID, true);
+ RemoteIntersection mIntersection;
+ mIntersection.set(nmBody, i, rIsID, true);
faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
rIsID++;
@@ -372,7 +392,6 @@ class SupportPatchFactory
void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, std::set<size_t>& patchBoundary, Patch& patchDofs) {
// insert local dofs
- const auto& lfe = cache_.get(fineElem.type());
const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID)->gridView();
const auto& indexSet = gridView.indexSet();
@@ -394,28 +413,35 @@ class SupportPatchFactory
addToPatch(isDofs, patchBoundary, patchDofs);
}
} else {
- size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
+ /* size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
if (fineIntersectionMapper_.count(faceIdx)) {
- const auto& intersections = coarseIntersectionMapper_[faceIdx];
+ const auto& intersections = fineIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
const auto& intersection = intersections[i];
- const auto& rIs = intersection.get();
+ const auto& glue = *fineContactNetwork_.glue(intersection.contactCouplingID);
+ const auto& rIs = intersection.get(glue);
const size_t outsideID = intersection.bodyID;
- std::set<size_t> rIsDofs;
+ std::vector<std::set<size_t>> rIsDofs;
remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs);
if (rIs.neighbor()) {
- auto outsideFather = coarseFather(rIs.outside(), fineContactNetwork_.body(outsideID)->gridView().level());
-
+ Element outsideFather;
+ if (intersection.flip)
+ outsideFather = coarseFather(rIs.inside(), fineContactNetwork_.body(outsideID)->level());
+ else
+ outsideFather = coarseFather(rIs.outside(), fineContactNetwork_.body(outsideID)->level());
+
+ std::set<size_t> totalRIsDofs(rIsDofs[0]);
+ totalRIsDofs.insert(rIsDofs[1].begin(), rIsDofs[1].end());
if (!coarsePatch.count(coarseIndices_.elementIndex(outsideID, outsideFather))) {
- addToPatchBoundary(rIsDofs, patchBoundary, patchDofs);
+ addToPatchBoundary(totalRIsDofs, patchBoundary, patchDofs);
} else {
- addToPatch(rIsDofs, patchBoundary, patchDofs);
+ addToPatch(totalRIsDofs, patchBoundary, patchDofs);
}
} else {
- addToPatchBoundary(rIsDofs, patchBoundary, patchDofs);
+ addToPatchBoundary(rIsDofs[0], patchBoundary, patchDofs);
}
}
} else {
@@ -423,12 +449,12 @@ class SupportPatchFactory
intersectionDofs(fineIndices_, is, bodyID, isDofs);
addToPatchBoundary(isDofs, patchBoundary, patchDofs);
- }
+ } */
}
}
}
- auto&& coarseFather(const Element& fineElem, const size_t coarseLevel) const {
+ auto coarseFather(const Element& fineElem, const size_t coarseLevel) const {
Element coarseElem = fineElem;
while (coarseElem.level() > coarseLevel) {
coarseElem = coarseElem.father();
@@ -443,7 +469,6 @@ class SupportPatchFactory
}
}
- template <class Intersection>
void addToPatch(const std::set<size_t>& dofs, const std::set<size_t>& patchBoundary, Patch& patchDofs) {
for (auto& dof : dofs) {
if (!patchBoundary.count(dof)) {
@@ -453,8 +478,13 @@ class SupportPatchFactory
}
template <class RIntersection>
- void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::array<std::set<size_t>, 2>& dofs) {
- dofs.clear();
+ void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::vector<std::set<size_t>>& dofs) {
+ dofs.resize(2);
+ dofs[0].clear();
+ dofs[1].clear();
+
+ const auto& isGeo = is.geometry();
+ const auto& isRefElement = Dune::ReferenceElements<ctype, dim-1>::general(is.type());
const auto& inside = is.inside();
const auto& insideGeo = inside.geometry();
@@ -468,18 +498,17 @@ class SupportPatchFactory
}
}
- const auto& isGeo = is.geometry();
- const auto& isRefElement = Dune::ReferenceElements<ctype, dim-1>::general(is.type());
-
- const auto& outside = is.outside();
- const auto& outsideGeo = outside.geometry();
- const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+ if (is.neighbor()) {
+ const auto& outside = is.outside();
+ const auto& outsideGeo = outside.geometry();
+ const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
- for (int i=0; i<outsideRefElement.size(dim); i++) {
- const auto& localCorner = isGeo.local(outsideGeo.corner(i));
+ for (int i=0; i<outsideRefElement.size(dim); i++) {
+ const auto& localCorner = isGeo.local(outsideGeo.corner(i));
- if (isRefElement.checkInside(localCorner)) {
- dofs[1].insert(indices.vertexIndex(outsideID, outside, i));
+ if (isRefElement.checkInside(localCorner)) {
+ dofs[1].insert(indices.vertexIndex(outsideID, outside, i));
+ }
}
}
}
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index 8ca00f64..b2004776 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
dune_add_test(SOURCES globalfrictioncontainertest.cc)
dune_add_test(SOURCES gridgluefrictiontest.cc)
dune_add_test(SOURCES nodalweightstest.cc)
+dune_add_test(SOURCES supportpatchfactorytest.cc)
diff --git a/src/tests/supportpatchfactorytest.cc b/src/tests/supportpatchfactorytest.cc
new file mode 100644
index 00000000..5d97ee91
--- /dev/null
+++ b/src/tests/supportpatchfactorytest.cc
@@ -0,0 +1,144 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#define MY_DIM 2
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/fufem/formatstring.hh>
+
+#include "../assemblers.hh"
+#include "../gridselector.hh"
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include "../factories/threeblocksfactory.cc"
+
+#include "../spatial-solving/preconditioners/supportpatchfactory.hh"
+#include "../utils/debugutils.hh"
+
+
+
+//#include <dune/tectonic/transformedglobalratestatefriction.hh>
+
+#include "common.hh"
+
+const int dim = MY_DIM;
+const int n = 5;
+const bool simplexGrid = true;
+
+const std::string path = "";
+const std::string outputFile = "supportpatchfactorytest.log";
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem-%dD.cfg", dim), parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+ return parset;
+}
+
+template <class SupportPatchFactory>
+void testSuite(const SupportPatchFactory& supportPatchFactory, const size_t bodyID, const int patchDepth = 0) {
+/* const auto& coarseContactNetwork = supportPatchFactory.coarseContactNetwork();
+ const auto& gridView = coarseContactNetwork.body(bodyID)->gridView();
+
+ for (const auto& e : elements(gridView)) {
+
+ }
+
+ using Patch = typename SupportPatchFactory::Patch;
+ Patch patch0, patch1, patch2, patch3;
+
+ // (const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0)
+
+ // interior patch inside of one body
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch0, patchDepth);
+
+ // patch at friction interface with two bodies in contact
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch1 patchDepth);
+
+ // patch at friction interface with two bodies in contact and with dirichlet boundary
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch2, patchDepth);
+
+ // crosspoint patch, where all 3 bodies are in contact
+ supportPatchFactory.build(0, const Element& coarseElement, const size_t localVertex, patch3, patchDepth);*/
+}
+
+int main(int argc, char *argv[]) { try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::ofstream out(path + outputFile);
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+ std::cout << "-------------------------------" << std::endl;
+ std::cout << "-- SupportPatchFactory Test: --" << std::endl;
+ std::cout << "-------------------------------" << std::endl << std::endl;
+
+ const auto parset = getParameters(argc, argv);
+
+ // set up contact network
+ ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
+ //using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
+ //threeBlocksFactory.build();
+
+ //ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork();
+ // threeBlocksFactory.build();
+
+ // auto& contactNetwork = threeBlocksFactory.contactNetwork();
+ //const size_t bodyCount = contactNetwork.nBodies();
+
+ //for (size_t i=0; i<bodyCount; i++) {
+ //printDofLocation(levelContactNetwork.leafView(i));
+
+ /*auto& levelViews = levelContactNetwork.levelViews(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"); */
+ //}
+
+ // assemble levelContactNetwork
+ //contactNetwork.assemble();
+
+ /* const auto& coarseContactNetwork = contactNetwork.level(0);
+ const auto& fineContactNetwork = contactNetwork.level(1);
+
+ SupportPatchFactory<decltype(coarseContactNetwork)> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
+
+
+ // set levelContactNetwork
+ LevelContactNetwork levelContactNetwork;
+ NBodyAssembler nBodyAssembler(grids.size(), 1);
+
+ std::vector<const GridType*> grids_ptr(grids.size());
+ for (size_t i=0; i<grids_ptr.size(); i++) {
+ grids_ptr[i] = grids[i].get();
+ }
+
+ nBodyAssembler.setGrids(grids_ptr);
+ nBodyAssembler.setCoupling(coupling, 0);
+
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();*/
+
+ bool passed = true;
+
+ std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+ 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;
+} // end try
+} // end main
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index 74396864..e2608f9b 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -17,8 +17,8 @@ void IterationRegister::reset() {
finalCount = FixedPointIterationCounter();
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeStepper(
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTimeStepper(
Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
const IgnoreVector& ignoreNodes,
@@ -28,7 +28,7 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeS
double relativeTime,
double relativeTau,
ExternalForces& externalForces,
- const std::vector<const ErrorNorm* >& errorNorms,
+ const ErrorNorms& errorNorms,
std::function<bool(Updaters &, Updaters &)> mustRefine)
: relativeTime_(relativeTime),
relativeTau_(relativeTau),
@@ -44,13 +44,13 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeS
mustRefine_(mustRefine),
errorNorms_(errorNorms) {}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-bool AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::reachedEnd() {
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+bool AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::reachedEnd() {
return relativeTime_ >= 1.0;
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::advance() {
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::advance() {
/*
| C | We check here if making the step R1 of size tau is a
| R1 | R2 | good idea. To check if we can coarsen, we compare
@@ -121,9 +121,9 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
return iterationRegister_;
}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-typename AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::UpdatersWithCount
-AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+typename AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(
Updaters const &oldUpdaters, double rTime, double rTau) {
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index f1f9b07c..af9fb0b1 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -15,7 +15,7 @@ struct IterationRegister {
FixedPointIterationCounter finalCount;
};
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
class AdaptiveTimeStepper {
struct UpdatersWithCount {
Updaters updaters;
@@ -27,7 +27,7 @@ class AdaptiveTimeStepper {
//using ConvexProblem = typename Factory::ConvexProblem;
//using Nonlinearity = typename Factory::Nonlinearity;
- using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>;
+ using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>;
using GlobalFriction = typename MyCoupledTimeStepper::GlobalFriction;
using BitVector = typename MyCoupledTimeStepper::BitVector;
@@ -44,7 +44,7 @@ class AdaptiveTimeStepper {
double relativeTime,
double relativeTau,
ExternalForces& externalForces,
- const std::vector<const ErrorNorm* >& errorNorms,
+ const ErrorNorms& errorNorms,
std::function<bool(Updaters &, Updaters &)> mustRefine);
bool reachedEnd();
@@ -69,7 +69,7 @@ class AdaptiveTimeStepper {
UpdatersWithCount R1_;
ExternalForces& externalForces_;
std::function<bool(Updaters &, Updaters &)> mustRefine_;
- const std::vector<const ErrorNorm* >& errorNorms_;
+ const ErrorNorms& errorNorms_;
IterationRegister iterationRegister_;
};
diff --git a/src/time-stepping/adaptivetimestepper_tmpl.cc b/src/time-stepping/adaptivetimestepper_tmpl.cc
index 7908b085..bda3c4c7 100644
--- a/src/time-stepping/adaptivetimestepper_tmpl.cc
+++ b/src/time-stepping/adaptivetimestepper_tmpl.cc
@@ -8,22 +8,22 @@
#include <dune/solvers/norms/energynorm.hh>
#include "../spatial-solving/solverfactory_tmpl.cc"
-#include "../data-structures/levelcontactnetwork_tmpl.cc"
+#include "../data-structures/contactnetwork_tmpl.cc"
#include "rate/rateupdater.hh"
#include "state/stateupdater.hh"
#include "updaters.hh"
-using NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
-using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
-using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
+using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
+using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
-template class AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
+template class AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index c5fbb89e..d3b9b1e2 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -4,15 +4,15 @@
#include "coupledtimestepper.hh"
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
-CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::CoupledTimeStepper(
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeStepper(
double finalTime, Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
Updaters updaters,
- const std::vector<const ErrorNorm* >& errorNorms,
+ const ErrorNorms& errorNorms,
ExternalForces& externalForces)
: finalTime_(finalTime),
parset_(parset),
@@ -24,9 +24,9 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::CoupledTimeSte
externalForces_(externalForces),
errorNorms_(errorNorms) {}
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
FixedPointIterationCounter
-CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(double relativeTime,
+CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double relativeTime,
double relativeTau) {
std::cout << "CoupledTimeStepper::step()" << std::endl;
@@ -48,7 +48,7 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(double re
updaters_.state_->setup(tau);
updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
- FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm> fixedPointIterator(
+ FixedPointIterator fixedPointIterator(
parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
velocityRHS, velocityIterate);
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index 4610c7a7..8b017775 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -8,16 +8,17 @@
#include "../spatial-solving/fixedpointiterator.hh"
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
class CoupledTimeStepper {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
using IgnoreVector = typename Factory::BitVector;
- //using Nonlinearity = typename Factory::Nonlinearity;
+ using FixedPointIterator = FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+
public:
- using GlobalFriction = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::GlobalFriction;
- using BitVector = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::BitVector;
- using ExternalForces = std::vector<const std::function<void(double, Vector &)>*>;
+ using GlobalFriction = typename FixedPointIterator::GlobalFriction;
+ using BitVector = typename FixedPointIterator::BitVector;
+ using ExternalForces = std::vector<std::unique_ptr<const std::function<void(double, Vector &)>>>;
public:
CoupledTimeStepper(double finalTime,
@@ -26,7 +27,8 @@ class CoupledTimeStepper {
const IgnoreVector& ignoreNodes,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
- Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
+ Updaters updaters,
+ const ErrorNorms& errorNorms,
ExternalForces& externalForces);
FixedPointIterationCounter step(double relativeTime, double relativeTau);
@@ -42,6 +44,6 @@ class CoupledTimeStepper {
Updaters updaters_;
ExternalForces& externalForces_;
- const std::vector<const ErrorNorm* >& errorNorms_;
+ const ErrorNorms& errorNorms_;
};
#endif
diff --git a/src/time-stepping/coupledtimestepper_tmpl.cc b/src/time-stepping/coupledtimestepper_tmpl.cc
index c70a229f..4d9f7f8e 100644
--- a/src/time-stepping/coupledtimestepper_tmpl.cc
+++ b/src/time-stepping/coupledtimestepper_tmpl.cc
@@ -8,21 +8,21 @@
#include <dune/solvers/norms/energynorm.hh>
#include "../spatial-solving/solverfactory_tmpl.cc"
-#include "../data-structures/levelcontactnetwork_tmpl.cc"
+#include "../data-structures/contactnetwork_tmpl.cc"
#include "rate/rateupdater.hh"
#include "state/stateupdater.hh"
#include "updaters.hh"
-using NBodyAssembler = typename MyLevelContactNetwork::NBodyAssembler;
-using BoundaryNodes = typename MyLevelContactNetwork::BoundaryNodes;
-using BoundaryFunctions = typename MyLevelContactNetwork::BoundaryFunctions;
+using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
-template class CoupledTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorm>;
+template class CoupledTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
diff --git a/src/time-stepping/rate/rateupdater_tmpl.cc b/src/time-stepping/rate/rateupdater_tmpl.cc
index aa8a7168..09299039 100644
--- a/src/time-stepping/rate/rateupdater_tmpl.cc
+++ b/src/time-stepping/rate/rateupdater_tmpl.cc
@@ -5,10 +5,10 @@
#include "../../explicitgrid.hh"
#include "../../explicitvectors.hh"
-#include "../../data-structures/levelcontactnetwork.hh"
+#include "../../data-structures/contactnetwork_tmpl.cc"
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
template class RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
template class Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index 9627a8ce..771a1a0c 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -1,8 +1,8 @@
-#include "../data-structures/levelcontactnetwork_tmpl.cc"
+#include "../data-structures/contactnetwork_tmpl.cc"
#include "../explicitvectors.hh"
-using BoundaryFunctions = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryFunctions;
-using BoundaryNodes = typename LevelContactNetwork<Grid, MY_DIM>::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
template
auto initRateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index 68ef4361..32c923f4 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -155,6 +155,20 @@
std::cout << message << " Step " << stepNumber << "!" << std::endl;
}*/
+ template <class GridView, class VectorType>
+ void writeToVTK(const GridView& gridView, const VectorType& x, const std::string path, const std::string name) {
+ Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+ std::ofstream lStream( "garbage.txt" );
+ std::streambuf* lBufferOld = std::cout.rdbuf();
+ std::cout.rdbuf(lStream.rdbuf());
+
+ vtkwriter.addVertexData(x,"data");
+ vtkwriter.pwrite(name, path, path);
+
+ std::cout.rdbuf( lBufferOld );
+ }
+
template <class GridView>
void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
Dune::VTKWriter<GridView> vtkwriter(gridView);
@@ -168,6 +182,7 @@
std::cout.rdbuf( lBufferOld );
}
+ /*
template <class VectorType, class DGBasisType>
void writeToVTK(const DGBasisType& basis, const VectorType& x, const std::string path, const std::string name) {
Dune::VTKWriter<typename DGBasisType::GridView> vtkwriter(basis.getGridView());
@@ -183,7 +198,7 @@
vtkwriter.pwrite(name, path, path);
std::cout.rdbuf( lBufferOld );
- }
+ }*/
template <class Intersection, class Basis>
void printIntersection(const Intersection& it, const Basis& basis, std::string message = "") {
--
GitLab