diff --git a/dune/tectonic/CMakeLists.txt b/dune/tectonic/CMakeLists.txt
index 7e6ac59121e6adcd6e5f28333526e99e862bf3c0..525c9d2bfe70b60425a80194ac6b5b808a9def07 100644
--- a/dune/tectonic/CMakeLists.txt
+++ b/dune/tectonic/CMakeLists.txt
@@ -1,5 +1,5 @@
install(FILES
- body.hh
+ bodydata.hh
frictiondata.hh
frictionpotential.hh
globalfrictiondata.hh
diff --git a/dune/tectonic/body.hh b/dune/tectonic/bodydata.hh
similarity index 85%
rename from dune/tectonic/body.hh
rename to dune/tectonic/bodydata.hh
index 2c2ed7214d673c3f91b1448463834634b111ffb9..eee19ff11ec6a919de1339ff09a29a0ae7238edd 100644
--- a/dune/tectonic/body.hh
+++ b/dune/tectonic/bodydata.hh
@@ -1,7 +1,7 @@
-#ifndef DUNE_TECTONIC_BODY_HH
-#define DUNE_TECTONIC_BODY_HH
+#ifndef DUNE_TECTONIC_BODY_DATA_HH
+#define DUNE_TECTONIC_BODY_DATA_HH
-template <int dimension> struct Body {
+template <int dimension> struct BodyData {
using ScalarFunction =
Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
Dune::FieldVector<double, 1>>;
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 3bb473174c2eff23faa233d711fc63793dca8b54..f2b048493b5f2de31b439f2023f7e211f8c76054 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,6 @@
set(SW_SOURCE_FILES
assemblers.cc
+ body.cc
enumparser.cc
hdf5/frictionalboundary-writer.cc
hdf5/iteration-writer.cc
@@ -22,7 +23,9 @@ set(SW_SOURCE_FILES
set(MSW_SOURCE_FILES
assemblers.cc
+ body.cc
enumparser.cc
+ levelcontactnetwork.cc
hdf5/frictionalboundary-writer.cc
hdf5/iteration-writer.cc
hdf5/patchinfo-writer.cc
@@ -34,6 +37,7 @@ set(MSW_SOURCE_FILES
multi-body-problem.cc
spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc
+ stackedblocksfactory.cc
time-stepping/adaptivetimestepper.cc
time-stepping/coupledtimestepper.cc
time-stepping/rate.cc
diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index 0d5ea0c609ac287f67b8df5e79277f17772ce894..036dbc9ee59e8e5c5d964e06d0477d693cbb71b1 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -3,5 +3,8 @@
#endif
#include "explicitgrid.hh"
+#include "explicitvectors.hh"
-template class MyAssembler<DeformedGrid::LeafGridView, MY_DIM>;
+#include "body.hh"
+
+template class MyAssembler<Body<Grid, Vector, MY_DIM>::DeformedLeafGridView, MY_DIM>;
diff --git a/src/body.cc b/src/body.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f110eec13b143b1b619a35c5ac4f248cef133785
--- /dev/null
+++ b/src/body.cc
@@ -0,0 +1,81 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+
+#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) :
+ 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_)),
+ matrices_()
+{
+ leafBoundaryConditions_.resize(0);
+
+ levelBoundaryConditions_.resize(grid_->maxLevel()+1);
+ for (size_t i=0; i<levelBoundaryConditions_.size(); i++) {
+ levelBoundaryConditions_[i].resize(0);
+ }
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assemble() {
+
+ // assemble matrices
+ assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
+ assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
+ assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
+
+ // assemble state energy norm
+ typename FrictionBoundaryCondition::BoundaryPatch frictionBoundary(this->leafView(), false); // boundary patch containing all friction boundaries (possibly disconnected)
+ for (size_t i=0; i<frictionBoundaryConditions_.size(); i++) {
+ const auto& frictionCondition = frictionBoundaryConditions_[i];
+ //std::cout << "grid Pointer " << frictionCondition->boundaryPatch()-> << std::endl;
+ frictionBoundary.addPatch(*frictionCondition->boundaryPatch());
+ }
+
+ ScalarMatrix relativeFrictionBoundaryMass;
+ assembler_->assembleFrictionalBoundaryMass(frictionBoundary, relativeFrictionBoundaryMass);
+ relativeFrictionBoundaryMass /= frictionBoundary.area(); // TODO: weight by individual friction patches?
+ stateEnergyNorm_ = std::make_shared<const EnergyNorm<ScalarMatrix, ScalarVector>>(relativeFrictionBoundaryMass);
+
+ // assemble forces
+ assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembleFriction(const Config::FrictionModel& frictionModel, const ScalarVector& weightedNormalStress) {
+ globalFriction_.resize(frictionBoundaryConditions_.size());
+
+ for (size_t i=0; i<globalFriction_.size(); i++) {
+ const auto& frictionBoundaryCondition = frictionBoundaryConditions_[i];
+ globalFriction_[i] = assembler_->assembleFrictionNonlinearity(frictionModel, *frictionBoundaryCondition->boundaryPatch(), *frictionBoundaryCondition->frictionData(), weightedNormalStress);
+ }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::externalForce(double relativeTime, VectorType& ell) const {
+ // Problem formulation: right-hand side
+ std::vector<std::shared_ptr<LeafBoundaryCondition>> leafNeumannConditions;
+ leafBoundaryConditions("neumann", leafNeumannConditions);
+ ell.resize(gravityFunctional_.size());
+ ell = gravityFunctional_;
+
+ for (size_t i=0; i<leafNeumannConditions.size(); i++) {
+ const auto& leafNeumannCondition = leafNeumannConditions[i];
+
+ VectorType b;
+ assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(), b, *leafNeumannCondition->boundaryFunction(),
+ relativeTime);
+ ell += b;
+ }
+}
+
+#include "body_tmpl.cc"
diff --git a/src/body.hh b/src/body.hh
new file mode 100644
index 0000000000000000000000000000000000000000..df86f2fa62f571ee793741666fac454609289870
--- /dev/null
+++ b/src/body.hh
@@ -0,0 +1,180 @@
+#ifndef SRC_BODY_HH
+#define SRC_BODY_HH
+
+#include <functional>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+//#include <dune/fufem/assemblers/assembler.hh>
+//#pragma clang diagnostic push
+//#pragma clang diagnostic ignored "-Wsign-compare"
+//#include <dune/fufem/functionspacebases/p0basis.hh>
+//#pragma clang diagnostic pop
+//#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/grid/geometrygrid/grid.hh>
+
+#include <dune/fufem/functions/deformationfunction.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+
+#include <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "assemblers.hh"
+#include "boundarycondition.hh"
+#include "enums.hh"
+#include "matrices.hh"
+
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims> class Body {
+public:
+ using GridType = GridTEMPLATE;
+ using VectorType = VectorTEMPLATE;
+
+ using DeformationFunction = DeformationFunction<typename GridType::LeafGridView, VectorType>;
+ using DeformedGridType = Dune::GeometryGrid<GridType, DeformationFunction>;
+ using DeformedLeafGridView = typename DeformedGridType::LeafGridView;
+ using DeformedLevelGridView = typename DeformedGridType::LevelGridView;
+
+ using Assembler = MyAssembler<DeformedLeafGridView, dims>;
+ using Matrix = typename Assembler::Matrix;
+ using LocalVector = typename VectorType::block_type;
+
+ using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+ using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+
+ using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
+ using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
+ using LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
+
+private:
+ std::shared_ptr<BodyData<dims>> bodyData_;
+
+ std::shared_ptr<GridType> grid_;
+ std::shared_ptr<DeformationFunction> deformation_;
+ std::shared_ptr<DeformedGridType> deformedGrid_;
+
+ DeformedLeafGridView leafView_;
+ int leafVertexCount_;
+
+ std::shared_ptr<Assembler> assembler_;
+ Matrices<Matrix,1> matrices_;
+
+ VectorType gravityFunctional_;
+ std::shared_ptr<const EnergyNorm<ScalarMatrix, ScalarVector>> stateEnergyNorm_;
+
+ // boundary conditions
+ std::vector<std::shared_ptr<LeafBoundaryCondition>> leafBoundaryConditions_;
+ std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>> levelBoundaryConditions_; // first index: level, second index: bc on lvl
+
+ // friction boundary conditions
+ std::vector<std::shared_ptr<FrictionBoundaryCondition>> frictionBoundaryConditions_;
+ std::vector<std::shared_ptr<GlobalFriction<Matrix, VectorType>>> globalFriction_;
+
+public:
+ Body(const std::shared_ptr<BodyData<dims>>& bodyData, const std::shared_ptr<GridType>& grid);
+
+ void assemble();
+ void assembleFriction(const Config::FrictionModel& frictionModel, const ScalarVector& weightedNormalStress);
+
+ const std::shared_ptr<Assembler>& assembler() const {
+ return assembler_;
+ }
+
+ const std::shared_ptr<BodyData<dims>>& data() const {
+ return bodyData_;
+ }
+
+ void addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition) {
+ leafBoundaryConditions_.push_back(leafBoundaryCondition);
+ }
+
+ const std::vector<std::shared_ptr<LeafBoundaryCondition>>& leafBoundaryConditions() const {
+ return leafBoundaryConditions_;
+ }
+
+ void leafBoundaryConditions(const std::string& tag, std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const {
+ selectedConditions.resize(0);
+ for (size_t i=0; i<leafBoundaryConditions_.size(); i++) {
+ if (leafBoundaryConditions_[i]->tag() == tag)
+ selectedConditions.push_back(leafBoundaryConditions_[i]);
+ }
+ }
+
+ void addBoundaryCondition(std::shared_ptr<FrictionBoundaryCondition> frictionBoundaryCondition) {
+ frictionBoundaryConditions_.push_back(frictionBoundaryCondition);
+ }
+
+ const std::vector<std::shared_ptr<FrictionBoundaryCondition>>& frictionBoundaryConditions() const {
+ return frictionBoundaryConditions_;
+ }
+
+ void addBoundaryCondition(std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition, int level) {
+ levelBoundaryConditions_[level].push_back(levelBoundaryCondition);
+ }
+
+ const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& levelBoundaryConditions() const {
+ return levelBoundaryConditions_;
+ }
+
+ void levelBoundaryConditions(const std::string& tag, std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const {
+ selectedConditions.resize(levelBoundaryConditions_.size());
+ for (size_t level=0; level<levelBoundaryConditions_.size(); level++) {
+ const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
+ for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
+ if (levelBoundaryConditions[i]->tag() == tag)
+ selectedConditions[level].push_back(levelBoundaryConditions[i]);
+ }
+ }
+ }
+
+ const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions(int level) const {
+ return levelBoundaryConditions_[level];
+ }
+
+ void levelBoundaryConditions(const std::string& tag, int level, std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const {
+ selectedConditions.resize(0);
+ const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
+ for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
+ if (levelBoundaryConditions[i]->tag() == tag)
+ selectedConditions.push_back(levelBoundaryConditions[i]);
+ }
+ }
+
+ const DeformedGridType* deformedGrid() const {
+ return deformedGrid_.get();
+ }
+
+ DeformationFunction& deformation() {
+ return *deformation_;
+ }
+
+ void setDeformation(const VectorType& def) {
+ deformation_->setDeformation(def);
+ }
+
+ const DeformedLeafGridView& leafView() const {
+ return leafView_;
+ }
+
+ int leafVertexCount() const {
+ return leafVertexCount_;
+ }
+
+ const Matrices<Matrix,1>& matrices() const {
+ return matrices_;
+ }
+
+ void externalForce(double relativeTime, VectorType& ell) const;
+};
+#endif
diff --git a/src/body_tmpl.cc b/src/body_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..25f890ffcff142320e4854264e6c2e9b564a832b
--- /dev/null
+++ b/src/body_tmpl.cc
@@ -0,0 +1,8 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "explicitgrid.hh"
+#include "explicitvectors.hh"
+
+template class Body<Grid, Vector, MY_DIM>;
diff --git a/src/boundarycondition.hh b/src/boundarycondition.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2ab5f2b758d1e86086f2275c7d21f9a3964e1ca9
--- /dev/null
+++ b/src/boundarycondition.hh
@@ -0,0 +1,116 @@
+#ifndef SRC_BOUNDARYCONDITION_HH
+#define SRC_BOUNDARYCONDITION_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+
+template <class GridView, int dims>
+class BoundaryCondition {
+public:
+ using BoundaryPatch = BoundaryPatch<GridView>;
+
+private:
+ using Function = Dune::VirtualFunction<double, double>;
+
+ const std::string tag_; // store type of boundary condition, e.g. dirichlet, neumann, friction, etc
+
+ std::shared_ptr<BoundaryPatch> boundaryPatch_;
+ std::shared_ptr<Dune::BitSetVector<dims>> boundaryNodes_;
+ std::shared_ptr<Function> boundaryFunction_;
+
+public:
+ BoundaryCondition(const std::string& tag = "") :
+ tag_(tag)
+ {}
+
+ BoundaryCondition(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function, const std::string& tag = "") :
+ tag_(tag),
+ boundaryPatch_(patch),
+ boundaryFunction_(function)
+ {}
+
+ void setBoundaryPatch(const GridView& gridView, std::shared_ptr<Dune::BitSetVector<dims>> nodes) {
+ boundaryNodes_ = nodes;
+ boundaryPatch_ = std::make_shared<BoundaryPatch>(gridView, *nodes);
+ }
+
+ void setBoundaryPatch(std::shared_ptr<BoundaryPatch> patch) {
+ boundaryPatch_ = patch;
+ }
+
+ void setBoundaryPatch(const BoundaryPatch& patch) {
+ boundaryPatch_ = std::make_shared<BoundaryPatch>(patch);
+ }
+
+ void setBoundaryFunction(std::shared_ptr<Function> function) {
+ boundaryFunction_ = function;
+ }
+
+ void set(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function) {
+ boundaryPatch_ = patch;
+ boundaryFunction_ = function;
+ }
+
+ const std::string& tag() const {
+ return tag_;
+ }
+
+ const std::shared_ptr<BoundaryPatch>& boundaryPatch() const {
+ return boundaryPatch_;
+ }
+
+ const std::shared_ptr<Dune::BitSetVector<dims>>& boundaryNodes() const {
+ return boundaryNodes_;
+ }
+
+ const std::shared_ptr<Function>& boundaryFunction() const {
+ return boundaryFunction_;
+ }
+};
+
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+template <class GridView, class LocalVectorType, int dims>
+class FrictionBoundaryCondition : public BoundaryCondition<GridView, dims>{
+private:
+ using Base = BoundaryCondition<GridView, dims>;
+ using LocalVector = LocalVectorType;
+
+ // friction data
+ std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
+ std::shared_ptr<GlobalFrictionData<dims>> frictionData_;
+
+public:
+ FrictionBoundaryCondition() :
+ Base("friction")
+ {}
+
+ FrictionBoundaryCondition(std::shared_ptr<typename Base::BoundaryPatch> patch, std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch, std::shared_ptr<GlobalFrictionData<dims>> frictionData) :
+ Base(patch, nullptr, "friction"),
+ weakeningPatch_(weakeningPatch),
+ frictionData_(frictionData)
+ {}
+
+ void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
+ weakeningPatch_ = weakeningPatch;
+ }
+
+ void setFrictionData(std::shared_ptr<GlobalFrictionData<dims>> frictionData) {
+ frictionData_ = frictionData;
+ }
+
+ const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch() const {
+ return weakeningPatch_;
+ }
+
+ const std::shared_ptr<GlobalFrictionData<dims>>& frictionData() const {
+ return frictionData_;
+ }
+};
+#endif
diff --git a/src/explicitgrid.hh b/src/explicitgrid.hh
index 685d2a49b434598a9be69ea4841e8251b55ab571..beefb64cc800f06604d3ac12484eb0be305e7902 100644
--- a/src/explicitgrid.hh
+++ b/src/explicitgrid.hh
@@ -11,4 +11,5 @@ using LevelGridView = Grid::LevelGridView;
using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
using DeformedGrid = DeformedGridComplex::DeformedGridType;
+
#endif
diff --git a/src/levelcontactnetwork.cc b/src/levelcontactnetwork.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8f58d4d250fa29abe3c268667db81a1fd9da675a
--- /dev/null
+++ b/src/levelcontactnetwork.cc
@@ -0,0 +1,125 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "levelcontactnetwork.hh"
+
+template <class GridType, int dimension>
+LevelContactNetwork<GridType, dimension>::LevelContactNetwork(int nBodies, int nCouplings, int level) :
+ level_(level),
+ bodies_(nBodies),
+ couplings_(nCouplings),
+ deformedGrids_(nBodies),
+ leafViews_(nBodies),
+ levelViews_(nBodies),
+ leafVertexCounts_(nBodies),
+ nBodyAssembler_(nBodies, nCouplings)
+{}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::assemble() {
+ for (size_t i=0; i<nBodies(); i++) {
+ bodies_[i]->assemble();
+ }
+
+ totalNodes("dirichlet", totalDirichletNodes_);
+
+ // set up dune-contact nBodyAssembler
+ nBodyAssembler_.setGrids(deformedGrids_);
+
+ for (size_t i=0; i<nCouplings(); i++) {
+ nBodyAssembler_.setCoupling(*couplings_[i], i);
+ }
+
+ nBodyAssembler_.assembleTransferOperator();
+ nBodyAssembler_.assembleObstacle();
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::assembleFriction(const Config::FrictionModel& frictionModel, const std::vector<ScalarVector>& weightedNormalStress) {
+ for (size_t i=0; i<nBodies(); i++) {
+ bodies_[i]->assembleFriction(frictionModel, weightedNormalStress[i]);
+ }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::setBodies(const std::vector<std::shared_ptr<Body>>& bodies) {
+ assert(nBodies()==bodies.size());
+ bodies_ = bodies;
+
+ matrices_.elasticity.resize(nBodies());
+ matrices_.damping.resize(nBodies());
+ matrices_.mass.resize(nBodies());
+
+ for (size_t i=0; i<nBodies(); i++) {
+ deformedGrids_[i] = bodies_[i]->deformedGrid();
+
+ leafViews_[i] = std::make_unique<DeformedLeafGridView>(deformedGrids_[i]->leafGridView());
+ levelViews_[i] = std::make_unique<DeformedLevelGridView>(deformedGrids_[i]->levelGridView(0));
+
+ leafVertexCounts_[i] = leafViews_[i]->size(dimension);
+
+ matrices_.elasticity[i] = bodies_[i]->matrices().elasticity;
+ matrices_.damping[i] = bodies_[i]->matrices().damping;
+ matrices_.mass[i] = bodies_[i]->matrices().mass;
+ }
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::setCouplings(const std::vector<std::shared_ptr<CouplingPair>>& couplings) {
+ assert(this->nCouplings()==couplings.size());
+ couplings_ = couplings;
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::constructBody(const std::shared_ptr<BodyData<dimension>>& bodyData,
+ const std::shared_ptr<GridType>& grid,
+ std::shared_ptr<Body>& body) const
+{
+ body = std::make_shared<Body>(bodyData, grid);
+}
+
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::constructCoupling(int nonMortarBodyIdx, int mortarBodyIdx,
+ const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
+ const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
+ std::shared_ptr<CouplingPair>& coupling) const
+{
+ coupling = std::make_shared<CouplingPair>();
+
+ auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<DeformedLeafBoundaryPatch>>();
+ std::shared_ptr<typename CouplingPair::BackEndType> backend = nullptr;
+ coupling->set(nonMortarBodyIdx, mortarBodyIdx, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
+}
+
+// collects all leaf boundary nodes from the different bodies identified by "tag" in a single BitSetVector
+template <class GridType, int dimension>
+void LevelContactNetwork<GridType, dimension>::totalNodes(const std::string& tag, Dune::BitSetVector<dimension>& totalNodes) {
+
+ int totalSize = 0;
+ for (size_t i=0; i<nBodies(); i++) {
+ totalSize += this->body(i)->leafVertexCount();
+ }
+
+ totalNodes.resize(totalSize);
+
+ int idx=0;
+ for (size_t i=0; i<nBodies(); i++) {
+ const auto& body = this->body(i);
+ std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+ body->leafBoundaryConditions(tag, boundaryConditions);
+
+ const int idxBackup = idx;
+ for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+ const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+ for (size_t j=0; j<nodes->size(); j++, idx++)
+ for (int k=0; k<dimension; k++)
+ totalNodes[idx][k] = (*nodes)[j][k];
+ idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+ }
+ }
+}
+
+#include "levelcontactnetwork_tmpl.cc"
diff --git a/src/levelcontactnetwork.hh b/src/levelcontactnetwork.hh
new file mode 100644
index 0000000000000000000000000000000000000000..be726c996f7cf0ca0300c4058534dbbf885047a2
--- /dev/null
+++ b/src/levelcontactnetwork.hh
@@ -0,0 +1,145 @@
+#ifndef SRC_LEVELCONTACTNETWORK_HH
+#define SRC_LEVELCONTACTNETWORK_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/bvector.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+//#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+
+#include "assemblers.hh"
+#include "body.hh"
+#include "enums.hh"
+#include "matrices.hh"
+//#include "multi-body-problem-data/myglobalfrictiondata.hh"
+
+template <class GridType, int dimension> class LevelContactNetwork {
+ private:
+ 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 CouplingPair = Dune::Contact::CouplingPair<DeformedGridType, DeformedGridType, field_type>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+
+ private:
+ const int level_;
+
+ std::vector<std::shared_ptr<Body>> bodies_;
+ std::vector<std::shared_ptr<CouplingPair>> couplings_;
+
+ std::vector<const DeformedGridType*> deformedGrids_;
+ std::vector<std::unique_ptr<const DeformedLeafGridView>> leafViews_;
+ std::vector<std::unique_ptr<const DeformedLevelGridView>> levelViews_;
+ std::vector<size_t> leafVertexCounts_;
+
+ Dune::Contact::NBodyAssembler<DeformedGridType, Vector> nBodyAssembler_;
+ Dune::BitSetVector<dimension> totalDirichletNodes_;
+
+ Matrices<Matrix,2> matrices_;
+ // std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionData_;
+ std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction_;
+
+ public:
+ LevelContactNetwork(int nBodies, int nCouplings, int level = 0);
+
+ void assemble();
+ void assembleFriction(const Config::FrictionModel& frictionModel, const std::vector<ScalarVector>& weightedNormalStress);
+
+ int level() const {
+ return level_;
+ }
+
+ void setBodies(const std::vector<std::shared_ptr<Body>>& bodies);
+ void setCouplings(const std::vector<std::shared_ptr<CouplingPair>>& couplings);
+
+ void constructBody(const std::shared_ptr<BodyData<dimension>>& bodyData,
+ const std::shared_ptr<GridType>& grid,
+ std::shared_ptr<Body>& body) const;
+ void constructCoupling(int nonMortarBodyIdx, int mortarBodyIdx,
+ const std::shared_ptr<DeformedLevelBoundaryPatch>& nonmortarPatch,
+ const std::shared_ptr<DeformedLevelBoundaryPatch>& mortarPatch,
+ std::shared_ptr<CouplingPair>& coupling) const;
+
+ void totalNodes(const std::string& tag, Dune::BitSetVector<dimension>& totalNodes);
+
+ size_t nBodies() const {
+ return bodies_.size();
+ }
+
+ size_t nCouplings() const {
+ return couplings_.size();
+ }
+
+ const std::vector<size_t>& leafVertexCounts() const {
+ return leafVertexCounts_;
+ }
+
+ const std::shared_ptr<Body>& body(int i) const {
+ return bodies_[i];
+ }
+
+ const std::vector<const DeformedGridType*>& deformedGrids() const {
+ return deformedGrids_;
+ }
+
+ const DeformedLeafGridView& leafView(int i) const {
+ return *leafViews_[i];
+ }
+
+ const DeformedLevelGridView& levelView(int i) const {
+ return *levelViews_[i];
+ }
+
+ int leafVertexCount(int i) const {
+ return leafVertexCounts_[i];
+ }
+
+ Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& nBodyAssembler() {
+ return nBodyAssembler_;
+ }
+
+ const Dune::Contact::NBodyAssembler<DeformedGridType, Vector>& nBodyAssembler() const {
+ return nBodyAssembler_;
+ }
+
+ const Matrices<Matrix,2>& matrices() const {
+ return matrices_;
+ }
+
+ const std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>>& globalFriction() const {
+ return globalFriction_;
+ }
+
+ const Dune::BitSetVector<dimension>& totalDirichletNodes() const {
+ return totalDirichletNodes_;
+ }
+};
+#endif
diff --git a/src/levelcontactnetwork_tmpl.cc b/src/levelcontactnetwork_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ee900a40719fa08df64fe9fb92838da6fc11fe05
--- /dev/null
+++ b/src/levelcontactnetwork_tmpl.cc
@@ -0,0 +1,7 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "explicitgrid.hh"
+
+template class LevelContactNetwork<Grid, MY_DIM>;
diff --git a/src/levelcontactnetworkfactory.hh b/src/levelcontactnetworkfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f9fb4e94ebd6c8396e12fe8ab49604150a1c5d86
--- /dev/null
+++ b/src/levelcontactnetworkfactory.hh
@@ -0,0 +1,54 @@
+#ifndef SRC_LEVELCONTACTNETWORKFACTORY_HH
+#define SRC_LEVELCONTACTNETWORKFACTORY_HH
+
+#include <dune/common/parametertree.hh>
+
+#include "levelcontactnetwork.hh"
+
+template <class GridType, int dims>
+class LevelContactNetworkFactory {
+public:
+ using LevelContactNetwork = LevelContactNetwork<GridType, dims>;
+
+protected:
+ using Body = typename LevelContactNetwork::Body;
+ using CouplingPair = typename LevelContactNetwork::CouplingPair;
+
+ const Dune::ParameterTree& parset_;
+ const size_t bodyCount_;
+ const size_t couplingCount_;
+
+ std::vector<std::shared_ptr<Body>> bodies_;
+ std::vector<std::shared_ptr<CouplingPair>> couplings_;
+
+ LevelContactNetwork levelContactNetwork_;
+
+private:
+ virtual void setBodies() = 0;
+ virtual void setCouplings() = 0;
+ virtual void setBoundaryConditions() = 0;
+
+public:
+ LevelContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
+ parset_(parset),
+ bodyCount_(bodyCount),
+ couplingCount_(couplingCount),
+ bodies_(bodyCount_),
+ couplings_(couplingCount_),
+ levelContactNetwork_(bodyCount_, couplingCount_) {}
+
+ void build() {
+ setBodies();
+ levelContactNetwork_.setBodies(bodies_);
+
+ setCouplings();
+ levelContactNetwork_.setCouplings(couplings_);
+
+ setBoundaryConditions();
+ }
+
+ LevelContactNetwork& levelContactNetwork() {
+ return levelContactNetwork_;
+ }
+};
+#endif
diff --git a/src/matrices.hh b/src/matrices.hh
index 01fc682dcbaff5b21b88a03f8bb958faa3c52be2..fa8b65d2681f15f549c13b336c27803bde90196f 100644
--- a/src/matrices.hh
+++ b/src/matrices.hh
@@ -1,13 +1,14 @@
#ifndef SRC_MATRICES_HH
#define SRC_MATRICES_HH
-template <class Matrix> class Matrices {
+template <class Matrix, size_t n>
+class Matrices {
public:
std::vector<std::shared_ptr<Matrix>> elasticity;
std::vector<std::shared_ptr<Matrix>> damping;
std::vector<std::shared_ptr<Matrix>> mass;
- Matrices(size_t n) {
+ Matrices() {
elasticity.resize(n);
damping.resize(n);
mass.resize(n);
@@ -19,4 +20,18 @@ template <class Matrix> class Matrices {
}
}
};
+
+template <class Matrix>
+class Matrices<Matrix, 1> {
+public:
+ std::shared_ptr<Matrix> elasticity;
+ std::shared_ptr<Matrix> damping;
+ std::shared_ptr<Matrix> mass;
+
+ Matrices() :
+ elasticity(std::make_shared<Matrix>()),
+ damping(std::make_shared<Matrix>()),
+ mass(std::make_shared<Matrix>())
+ {}
+};
#endif
diff --git a/src/multi-body-problem-data/bc.hh b/src/multi-body-problem-data/bc.hh
index 7c29cf087571f8aca3e98d987ad3e711e991122d..ebc6f744a34a261bd4b829a47679b294e3b36cb8 100644
--- a/src/multi-body-problem-data/bc.hh
+++ b/src/multi-body-problem-data/bc.hh
@@ -1,6 +1,8 @@
#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#include <dune/common/function.hh>
+
class VelocityDirichletCondition
: public Dune::VirtualFunction<double, double> {
void evaluate(double const &relativeTime, double &y) const {
diff --git a/src/multi-body-problem-data/cuboidgeometry.cc b/src/multi-body-problem-data/cuboidgeometry.cc
index aaf4dfe0683a53ca567f50affd2840afcb1cdd6d..77cd62f9e7636f5c7bed8e2b150973cbf92c9f50 100644
--- a/src/multi-body-problem-data/cuboidgeometry.cc
+++ b/src/multi-body-problem-data/cuboidgeometry.cc
@@ -35,26 +35,36 @@ void CuboidGeometry::setupWeak(const LocalVector& weakOrigin, const double weakL
*/
#if MY_DIM == 3
-CuboidGeometry::CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength, const double length, const double width, const double depth_):
+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):
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}),
- X(weakOrigin),
- Y({X[0]+weakLength, X[1], 0}),
+ V(lowerWeakOrigin),
+ W({V[0]+lowerWeakLength*lengthScale, V[1], 0}),
+ X(upperWeakOrigin),
+ Y({X[0]+upperWeakLength*lengthScale, X[1], 0}),
depth(depth_*lengthScale) {}
#else
-CuboidGeometry::CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength, const double length, const double width):
+CuboidGeometry::CuboidGeometry(const LocalVector& origin,
+ const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
+ const double lowerWeakLength, const double upperWeakLength,
+ const double length, const double width):
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_}),
- X(weakOrigin),
- Y({X[0]+weakLength, X[1]}) {}
+ V(lowerWeakOrigin),
+ W({V[0]+lowerWeakLength*lengthScale, V[1]}),
+ X(upperWeakOrigin),
+ Y({X[0]+upperWeakLength*lengthScale, X[1]}) {}
#endif
@@ -64,10 +74,13 @@ void CuboidGeometry::write() const {
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;
}
+/*
void CuboidGeometry::render() const {
#ifdef HAVE_CAIROMM
std::string const filename = "geometry.png";
@@ -188,3 +201,4 @@ void CuboidGeometry::render() const {
surface->write_to_png(filename);
#endif
}
+*/
diff --git a/src/multi-body-problem-data/cuboidgeometry.hh b/src/multi-body-problem-data/cuboidgeometry.hh
index 242454483d07285ad0d2a32332884dbee9572b34..d189502eda9e007c7ebab67697c2e8a3eaba4eb7 100644
--- a/src/multi-body-problem-data/cuboidgeometry.hh
+++ b/src/multi-body-problem-data/cuboidgeometry.hh
@@ -14,64 +14,40 @@ class CuboidGeometry {
private:
const double length_;
const double width_;
- //const double weakLength_; // length of the weak zone
+ // const double lowerWeakLength_; // length of the lower weak zone
+ // const double upperWeakLength_; // length of the lower weak zone
public:
+ // corners of cube
const LocalVector A;
const LocalVector B;
const LocalVector C;
const LocalVector D;
+
+ // corners of lower weakening patch
+ const LocalVector V;
+ const LocalVector W;
+
+ // corners of upper weakening patch
const LocalVector X;
const LocalVector Y;
#if MY_DIM == 3
const double depth;
- CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength = 0.20, const double length = 1.00, const double width = 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);
#else
- CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength = 0.20, const double length = 1.00, const double width = 0.27);
+ 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);
#endif
void write() const;
- void render() const;
+ // void render() const;
};
#endif
-
- /* from Elias
- double const s = 1.0; // scaling factor
-
- double const rightLeg = 0.27 * s;
- double const leftLeg = 1.00 * s;
- double const leftAngle = atan(rightLeg / leftLeg);
- double const viscoHeight = 0.06 * s; // Height of the viscous bottom layer
- double const weakLen = 0.20 * s; // Length of the weak zone
-
- double const zDistance = 0.35;
-
- LocalVector2D const A = {0, 0};
- LocalVector2D const B = {leftLeg, -rightLeg};
- LocalVector2D const C = {leftLeg, 0};
-
- LocalVector2D const Z = {zDistance * s, 0};
- LocalVector2D const Y = {zDistance * s, -zDistance *s / leftLeg *rightLeg};
- LocalVector2D const X = {Y[0] - weakLen * std::cos(leftAngle),
- Y[1] + weakLen *std::sin(leftAngle)};
-
- LocalVector2D const U = {X[0], 0};
-
- LocalVector2D const K = {B[0] - leftLeg * viscoHeight / rightLeg,
- B[1] + viscoHeight};
- LocalVector2D const M = {B[0], B[1] + viscoHeight};
-
- LocalVector2D const G = midPoint(A, X);
- LocalVector2D const H = midPoint(X, Y);
- LocalVector2D const J = midPoint(Y, B);
-
- LocalVector2D const I = {Y[0] + G[0], Y[1] + G[1]};
-
- LocalVector2D const zenith = {0, 1};
-
- LocalMatrix2D const rotation = {{std::cos(leftAngle), -std::sin(leftAngle)},
- {std::sin(leftAngle), std::cos(leftAngle)}};
- */
diff --git a/src/multi-body-problem-data/mybody.hh b/src/multi-body-problem-data/mybody.hh
index 08b5d4b9bd9d177fdb82965c7ac2de7ac7b05d3e..a79e359ee7e6670fed506bd4a64c2885cad7e5e8 100644
--- a/src/multi-body-problem-data/mybody.hh
+++ b/src/multi-body-problem-data/mybody.hh
@@ -1,25 +1,26 @@
-#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODY_HH
-#define SRC_MULTI_BODY_PROBLEM_DATA_MYBODY_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
#include <dune/common/fvector.hh>
#include <dune/fufem/functions/constantfunction.hh>
-#include <dune/tectonic/body.hh>
+#include <dune/tectonic/bodydata.hh>
#include <dune/tectonic/gravity.hh>
#include "cuboidgeometry.hh"
#include "segmented-function.hh"
-template <int dimension> class MyBody : public Body<dimension> {
- using typename Body<dimension>::ScalarFunction;
- using typename Body<dimension>::VectorField;
+template <int dimension> class MyBodyData : public BodyData<dimension> {
+ using typename BodyData<dimension>::ScalarFunction;
+ using typename BodyData<dimension>::VectorField;
public:
- MyBody(Dune::ParameterTree const &parset)
+ MyBodyData(Dune::ParameterTree const &parset, const Dune::FieldVector<double, dimension>& zenith)
: poissonRatio_(parset.get<double>("body.poissonRatio")),
youngModulus_(3.0 * parset.get<double>("body.bulkModulus") *
(1.0 - 2.0 * poissonRatio_)),
+ zenith_(zenith),
shearViscosityField_(
parset.get<double>("body.elastic.shearViscosity"),
parset.get<double>("body.viscoelastic.shearViscosity")),
@@ -28,7 +29,7 @@ template <int dimension> class MyBody : public Body<dimension> {
parset.get<double>("body.viscoelastic.bulkViscosity")),
densityField_(parset.get<double>("body.elastic.density"),
parset.get<double>("body.viscoelastic.density")),
- gravityField_(densityField_, MyGeometry::zenith,
+ gravityField_(densityField_, zenith_,
parset.get<double>("gravity")) {}
double getPoissonRatio() const override { return poissonRatio_; }
@@ -47,6 +48,8 @@ template <int dimension> class MyBody : public Body<dimension> {
private:
double const poissonRatio_;
double const youngModulus_;
+ Dune::FieldVector<double, dimension> zenith_;
+
SegmentedFunction const shearViscosityField_;
SegmentedFunction const bulkViscosityField_;
SegmentedFunction const densityField_;
diff --git a/src/multi-body-problem-data/myglobalfrictiondata.hh b/src/multi-body-problem-data/myglobalfrictiondata.hh
index d92e7223b20d6e257b5dc3168169b93cdbfa8eb2..7666a467ccfff0472f6bd71d5fe76f27d7441800 100644
--- a/src/multi-body-problem-data/myglobalfrictiondata.hh
+++ b/src/multi-body-problem-data/myglobalfrictiondata.hh
@@ -1,5 +1,5 @@
-#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
#include <dune/common/function.hh>
@@ -14,14 +14,15 @@ class MyGlobalFrictionData : public GlobalFrictionData<LocalVector::dimension> {
public:
MyGlobalFrictionData(Dune::ParameterTree const &parset,
- ConvexPolyhedron<LocalVector> const &segment)
+ ConvexPolyhedron<LocalVector> const &segment,
+ double lengthScale = 1.0)
: C_(parset.get<double>("C")),
L_(parset.get<double>("L")),
V0_(parset.get<double>("V0")),
a_(parset.get<double>("strengthening.a"),
- parset.get<double>("weakening.a"), segment),
+ parset.get<double>("weakening.a"), segment, lengthScale),
b_(parset.get<double>("strengthening.b"),
- parset.get<double>("weakening.b"), segment),
+ parset.get<double>("weakening.b"), segment, lengthScale),
mu0_(parset.get<double>("mu0")) {}
double const &C() const override { return C_; }
diff --git a/src/multi-body-problem-data/patchfunction.hh b/src/multi-body-problem-data/patchfunction.hh
index 72cdc867cffd4eba8091a0a2da0fa366e196627d..1f7e9b08f007ad498025922417cd1cf93f47a93c 100644
--- a/src/multi-body-problem-data/patchfunction.hh
+++ b/src/multi-body-problem-data/patchfunction.hh
@@ -1,5 +1,5 @@
-#ifndef SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
@@ -17,13 +17,15 @@ class PatchFunction
double const v2_;
Polyhedron const &segment_;
+ double const lengthScale_;
+
public:
- PatchFunction(double v1, double v2, Polyhedron const &segment)
- : v1_(v1), v2_(v2), segment_(segment) {}
+ PatchFunction(double v1, double v2, Polyhedron const &segment, double lengthScale = 1.0)
+ : v1_(v1), v2_(v2), segment_(segment), lengthScale_(lengthScale) {}
void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
Dune::FieldVector<double, 1> &y) const {
- y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+ y = distance(x, segment_, 1e-6 * lengthScale_) <= 1e-5 ? v2_
: v1_;
}
};
diff --git a/src/multi-body-problem-data/segmented-function.hh b/src/multi-body-problem-data/segmented-function.hh
index cab2dfef6e7bd5b764bd04048e9f43d9593b2ba9..a75a7199559a5a47c8b8527a9a9fef13a6e842f0 100644
--- a/src/multi-body-problem-data/segmented-function.hh
+++ b/src/multi-body-problem-data/segmented-function.hh
@@ -17,7 +17,7 @@ class SegmentedFunction
return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
};
bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
- return liesBelow(MyGeometry::K, MyGeometry::M, z);
+ return false; //TODO liesBelow(MyGeometry::K, MyGeometry::M, z);
};
double const _v1;
diff --git a/src/multi-body-problem-data/weakpatch.hh b/src/multi-body-problem-data/weakpatch.hh
index e48f5725a93ac97f8e0fbfd1bfde0532e3f292d6..4e54105bc293a3142033c1e34b04f58ad58c9540 100644
--- a/src/multi-body-problem-data/weakpatch.hh
+++ b/src/multi-body-problem-data/weakpatch.hh
@@ -1,34 +1,67 @@
#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 "cuboidgeometry.hh"
template <class LocalVector>
-ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry) {
- ConvexPolyhedron<LocalVector> weakPatch;
+void getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry, ConvexPolyhedron<LocalVector>& lowerWeakPatch, ConvexPolyhedron<LocalVector>& upperWeakPatch) {
+
#if MY_DIM == 3
- weakPatch.vertices.resize(4);
- weakPatch.vertices[0] = weakPatch.vertices[2] = cuboidGeometry.X;
- weakPatch.vertices[1] = weakPatch.vertices[3] = cuboidGeometry.Y;
- for (size_t k = 0; k < 2; ++k) {
- weakPatch.vertices[k][2] = -cuboidGeometry.depth / 2.0;
- weakPatch.vertices[k + 2][2] = cuboidGeometry.depth / 2.0;
+ 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);
+ }
}
- switch (parset.get<Config::PatchType>("patchType")) {
- case Config::Rectangular:
- break;
- case Config::Trapezoidal:
- weakPatch.vertices[1][0] += 0.05 * cuboidGeometry.lengthScale;
- weakPatch.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
- weakPatch.vertices.resize(2);
- weakPatch.vertices[0] = cuboidGeometry.X;
- weakPatch.vertices[1] = cuboidGeometry.Y;
+ 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
- return weakPatch;
};
#endif
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index a352da55283f59f131df1d88a7c04bf16e820273..ca941176059bcbe7a0280ebcb286639b829ea868 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -57,15 +57,16 @@
#include "gridselector.hh"
#include "hdf5-writer.hh"
#include "hdf5/restart-io.hh"
+#include "levelcontactnetwork.hh"
#include "matrices.hh"
#include "program_state.hh"
#include "multi-body-problem-data/bc.hh"
#include "multi-body-problem-data/mybody.hh"
-#include "multi-body-problem-data/cuboidgeometry.hh"
-#include "multi-body-problem-data/myglobalfrictiondata.hh"
+
+
#include "multi-body-problem-data/mygrids.hh"
-#include "multi-body-problem-data/weakpatch.hh"
//#include "spatial-solving/solverfactory.hh"
+#include "stackedblocksfactory.hh"
//#include "time-stepping/adaptivetimestepper.hh"
#include "time-stepping/rate.hh"
#include "time-stepping/state.hh"
@@ -81,6 +82,20 @@
size_t const dims = MY_DIM;
+template <class GridView>
+void writeToVTK(const GridView& gridView, 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.pwrite(name, path, path);
+
+ std::cout.rdbuf( lBufferOld );
+}
+
+
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);
@@ -120,257 +135,35 @@ int main(int argc, char *argv[]) {
using ScalarVector = Assembler::ScalarVector;
using field_type = Matrix::field_type;
- // set up material properties of all bodies
- MyBody<dims> const body(parset);
-
- // set up cuboid geometries
- const size_t bodyCount = parset.get<int>("problem.bodyCount");
- std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries(bodyCount);
-
- double const length = 1.00;
- double const width = 0.27;
- double const weakLength = 0.20;
-
-#if MY_DIM == 3
- double const depth = 0.60;
- // TODO: replace with make_shared
- cuboidGeometries[0] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry({0,0,0}, {0.2,width,0}, weakLength, length, width, depth));
- for (size_t i=1; i<bodyCount; i++) {
- cuboidGeometries[i] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry(cuboidGeometries[i-1]->D, {0.6,i*width,0}, weakLength, length, width, depth));
- }
-#else
- // TODO: replace with make_shared
- cuboidGeometries[0] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry({0,0}, {0.2,width}, weakLength, length, width));
- for (size_t i=1; i<bodyCount; i++) {
- cuboidGeometries[i] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry(cuboidGeometries[i-1]->D, {0.6,i*width}, weakLength, length, width));
- }
-#endif
-
- // ------------
- // set up grids
- // ------------
- GridsConstructor<Grid> gridConstructor(cuboidGeometries);
- auto& grids = gridConstructor.getGrids();
-
- std::vector<ConvexPolyhedron<LocalVector>> weakPatches(grids.size());
-
- for (size_t i=0; i<grids.size(); i++) {
- const auto& cuboidGeometry = *cuboidGeometries[i];
-
- // define weak patch and refine grid
- auto& weakPatch = weakPatches[i];
- weakPatch = getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"), cuboidGeometry);
- refine(*grids[i], weakPatch, parset.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;
- }
-
- // set up deformed grids
- DeformedGridComplex deformedGridComplex;
- for (size_t i=0; i<grids.size(); i++) {
- deformedGridComplex.addGrid(*grids[i], "body" + std::to_string(i));
- }
-
- const auto deformedGrids = deformedGridComplex.gridVector();
-
- // ------------------------
- // set up contact couplings
- // ------------------------
- std::vector<std::shared_ptr<DefLeafGridView>> leafViews(deformedGrids.size());
- std::vector<std::shared_ptr<DefLevelGridView>> levelViews(deformedGrids.size());
- std::vector<size_t> leafVertexCounts(deformedGrids.size());
-
- using LeafFaces = MyFaces<DefLeafGridView>;
- using LevelFaces = MyFaces<DefLevelGridView>;
- std::vector<std::shared_ptr<LeafFaces>> leafFaces(deformedGrids.size());
- std::vector<std::shared_ptr<LevelFaces>> levelFaces(deformedGrids.size());
-
- for (size_t i=0; i<deformedGrids.size(); i++) {
- leafViews[i] = std::make_shared<DefLeafGridView>(deformedGrids[i]->leafGridView());
- levelViews[i] = std::make_shared<DefLevelGridView>(deformedGrids[i]->levelGridView(0));
-
- leafVertexCounts[i] = leafViews[i]->size(dims);
-
- const auto& cuboidGeometry = *cuboidGeometries[i];
- leafFaces[i] = std::make_shared<LeafFaces>(*leafViews[i], cuboidGeometry);
- levelFaces[i] = std::make_shared<LevelFaces>(*levelViews[i], cuboidGeometry);
- }
-
- // set up contact couplings
- using LeafBoundaryPatch = BoundaryPatch<DefLeafGridView>;
- using LevelBoundaryPatch = BoundaryPatch<DefLevelGridView>;
-
- using CouplingPair = Dune::Contact::CouplingPair<DeformedGrid, DeformedGrid, field_type>;
- std::vector<std::shared_ptr<CouplingPair>> couplings(bodyCount-1);
-
- for (size_t i=0; i<couplings.size(); i++) {
- couplings[i] = std::make_shared<CouplingPair>();
-
- auto nonmortarPatch = std::make_shared<LevelBoundaryPatch>(levelFaces[i]->upper);
- auto mortarPatch = std::make_shared<LevelBoundaryPatch>(levelFaces[i+1]->lower);
-
- auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
- std::shared_ptr<CouplingPair::BackEndType> backend = nullptr;
- couplings[i]->set(i, i+1, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
- }
-
- // ----------------------------------
- // set up dune-contact nBodyAssembler
- // ----------------------------------
- Dune::Contact::NBodyAssembler<DeformedGrid, Vector> nBodyAssembler(bodyCount, bodyCount-1);
- nBodyAssembler.setGrids(deformedGrids);
-
- for (size_t i=0; i<couplings.size(); i++) {
- nBodyAssembler.setCoupling(*couplings[i], i);
- }
-
- nBodyAssembler.assembleTransferOperator();
- nBodyAssembler.assembleObstacle();
-
- // --------------------------
- // set up boundary conditions
- // --------------------------
- std::vector<BoundaryPatch<DefLeafGridView>> neumannBoundaries(bodyCount);
- std::vector<BoundaryPatch<DefLeafGridView>> surfaces(bodyCount);
-
- // friction boundary
- std::vector<BoundaryPatch<DefLeafGridView>> frictionalBoundaries(bodyCount);
- std::vector<Dune::BitSetVector<1>> frictionalNodes(bodyCount);
-
- // Dirichlet boundary
- std::vector<Dune::BitSetVector<dims>> noNodes(bodyCount);
- std::vector<std::vector<Dune::BitSetVector<dims>>> dirichletNodes(bodyCount);
-
- // set up functions for time-dependent boundary conditions
- using Function = Dune::VirtualFunction<double, double>;
- std::vector<const Function*> velocityDirichletFunctions(grids.size());
- const Function& neumannFunction = NeumannCondition();
-
- for (size_t i=0; i<grids.size(); i++) {
- const auto& leafVertexCount = leafVertexCounts[i];
-
- std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
-
- // Neumann boundary
- neumannBoundaries[i] = BoundaryPatch<DefLeafGridView>(*leafViews[i]);
-
- // friction boundary
- auto& frictionalBoundary = frictionalBoundaries[i];
- if (i==0) {
- frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->upper);
- } else if (i==bodyCount-1) {
- frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->lower);
- } else {
- frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->lower);
- frictionalBoundary.addPatch(BoundaryPatch<DefLeafGridView>(leafFaces[i]->upper));
- }
- frictionalNodes[i] = *frictionalBoundary.getVertices();
- //surfaces[i] = BoundaryPatch<GridView>(myFaces.upper);
+ // ----------------------
+ // set up contact network
+ // ----------------------
+ StackedBlocksFactory<Grid, dims> stackedBlocksFactory(parset);
+ using LevelContactNetwork = typename StackedBlocksFactory<Grid, dims>::LevelContactNetwork;
+ stackedBlocksFactory.build();
- // TODO: Dirichlet Boundary
- velocityDirichletFunctions[i] = new VelocityDirichletCondition();
+ LevelContactNetwork& levelContactNetwork = stackedBlocksFactory.levelContactNetwork();
+ const size_t bodyCount = levelContactNetwork.nBodies();
- noNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
+ for (size_t i=0; i<levelContactNetwork.deformedGrids().size(); i++) {
+ /* Vector def(levelContactNetwork.deformedGrids()[i]->size(dims));
+ def = 1;
+ deformedGridComplex.setDeformation(def, i);*/
- auto& gridDirichletNodes = dirichletNodes[i];
- gridDirichletNodes = Dune::BitSetVector<dims>(leafVertexCount);
-
- for (int j=0; j<leafVertexCount; j++) {
- if (leafFaces[i]->right.containsVertex(j))
- gridDirichletNodes[j][0] = true;
-
- if (leafFaces[i]->lower.containsVertex(j))
- gridDirichletNodes[j][0] = true;
-
- #if MY_DIM == 3
- if (leafFaces[i]->front.containsVertex(j) || leafFaces[i]->back.containsVertex(j))
- gridDirichletNodes[j][2] = true;
- #endif
- }
+ writeToVTK(levelContactNetwork.leafView(i), "", "deformedGrid" + std::to_string(i));
}
- // extend Dirichlet information to all grid levels
- /* BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[i]);
- dirichletNodes[i].resize(toplevel+1);
-
- for (int j=0; j<=toplevel; j++) {
-
- int fSSize = grids[i]->size(j,dim);
- dirichletNodes[i][j].resize(fSSize);
- for (int k=0; k<fSSize; k++)
- for (int l=0; l<dim; l++)
- dirichletNodes[i][j][k][l] = dirichletBoundary[i][j].containsVertex(k);
-
- }
-
- sampleOnBitField(*grids[i], dirichletValues[i], dirichletNodes[i]);*/
-
- const int toplevel = deformedGrids[0]->maxLevel();
- std::vector<Dune::BitSetVector<dims> > totalDirichletNodes(toplevel+1);
+ // ----------------------------
+ // assemble levelContactNetwork
+ // ----------------------------
- for (int i=0; i<=toplevel; i++) {
- int totalSize = 0;
- for (int j=0; j<deformedGrids.size(); j++)
- totalSize += deformedGrids[j]->size(i, dims);
-
- totalDirichletNodes[i].resize(totalSize);
-
- int idx=0;
- for (int j=0; j<deformedGrids.size(); j++)
- for (size_t k=0; k<dirichletNodes[j][i].size(); k++)
- for (int l=0; l<dims; l++)
- totalDirichletNodes[i][idx++][l] = dirichletNodes[j][i][k][l];
- }
-
- // ------------------------------------------------------------------------------------
- // set up individual assemblers for each body, assemble problem (matrices, forces, rhs)
- // ------------------------------------------------------------------------------------
- std::vector<std::shared_ptr<Assembler>> assemblers(bodyCount);
- Matrices<Matrix> matrices(bodyCount);
-
- std::vector<Vector> gravityFunctionals(bodyCount);
- std::vector<std::function<void(double, Vector&)>> externalForces(bodyCount);
- std::vector<const EnergyNorm<ScalarMatrix, ScalarVector>* > stateEnergyNorms(bodyCount);
-
- for (size_t i=0; i<assemblers.size(); i++) {
- auto& assembler = assemblers[i];
- assembler = std::make_shared<Assembler>(*leafViews[i]);
-
- assembler->assembleElasticity(body.getYoungModulus(), body.getPoissonRatio(), *matrices.elasticity[i]);
- assembler->assembleViscosity(body.getShearViscosityField(), body.getBulkViscosityField(), *matrices.damping[i]);
- assembler->assembleMass(body.getDensityField(), *matrices.mass[i]);
-
- ScalarMatrix relativeFrictionalBoundaryMass;
- assembler->assembleFrictionalBoundaryMass(frictionalBoundaries[i], relativeFrictionalBoundaryMass);
- relativeFrictionalBoundaryMass /= frictionalBoundaries[i].area();
- stateEnergyNorms[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(relativeFrictionalBoundaryMass);
-
- // assemble forces
- assembler->assembleBodyForce(body.getGravityField(), gravityFunctionals[i]);
-
- // Problem formulation: right-hand side
- externalForces[i] =
- [&](double _relativeTime, Vector &_ell) {
- assemblers[i]->assembleNeumann(neumannBoundaries[i], _ell, neumannFunction,
- _relativeTime);
- _ell += gravityFunctionals[i];
- };
- }
+ levelContactNetwork.assemble();
// -----------------
// init input/output
// -----------------
+ const auto& leafVertexCounts = levelContactNetwork.leafVertexCounts();
using MyProgramState = ProgramState<Vector, ScalarVector>;
MyProgramState programState(leafVertexCounts);
auto const firstRestart = parset.get<size_t>("io.restarts.first");
@@ -399,36 +192,31 @@ int main(int argc, char *argv[]) {
if (firstRestart > 0) // automatically adjusts the time and timestep
restartIO->read(firstRestart, programState);
else
- programState.setupInitialConditions(parset, nBodyAssembler, externalForces,
- matrices, assemblers, totalDirichletNodes,
- noNodes, frictionalBoundaries, body);
-
-
- // assemble friction
- std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionInfo(weakPatches.size());
- std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction(weakPatches.size());
- for (size_t i=0; i<frictionInfo.size(); i++) {
- frictionInfo[i] = std::make_shared<MyGlobalFrictionData<LocalVector>>(parset.sub("boundary.friction"), weakPatches[i]);
- globalFriction[i] = assemblers[i]->assembleFrictionNonlinearity(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundaries[i], *frictionInfo[i], programState.weightedNormalStress[i]);
- globalFriction[i]->updateAlpha(programState.alpha[i]);
- }
+ programState.setupInitialConditions(parset, levelContactNetwork);
+
+ // ------------------------
+ // assemble global friction
+ // ------------------------
+ levelContactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+ /*
using MyVertexBasis = typename Assembler::VertexBasis;
using MyCellBasis = typename Assembler::CellBasis;
- std::vector<Vector> vertexCoordinates(leafVertexCounts.size());
- std::vector<const MyVertexBasis* > vertexBases(leafVertexCounts.size());
- std::vector<const MyCellBasis* > cellBases(leafVertexCounts.size());
+ std::vector<Vector> vertexCoordinates(bodyCount);
+ std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+ std::vector<const MyCellBasis* > cellBases(bodyCount);
- for (size_t i=0; i<vertexCoordinates.size(); i++) {
- vertexBases[i] = &(assemblers[i]->vertexBasis);
- cellBases[i] = &(assemblers[i]->cellBasis);
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = levelContactNetwork.body(i);
+ vertexBases[i] = &(body->assembler()->vertexBasis);
+ cellBases[i] = &(body->assembler()->cellBasis);
auto& vertexCoords = vertexCoordinates[i];
vertexCoords.resize(leafVertexCounts[i]);
Dune::MultipleCodimMultipleGeomTypeMapper<
- DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(*leafViews[i], Dune::mcmgVertexLayout());
- for (auto &&v : vertices(*leafViews[i]))
+ DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->leafView(), Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(body->leafView()))
vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
@@ -446,7 +234,7 @@ int main(int argc, char *argv[]) {
auto const report = [&](bool initial = false) {
if (writeData) {
- dataWriter->reportSolution(programState, globalFriction);
+ dataWriter->reportSolution(programState, levelContactNetwork.globalFriction());
if (!initial)
dataWriter->reportIterations(programState, iterationCount);
dataFile->flush();
@@ -465,11 +253,12 @@ int main(int argc, char *argv[]) {
<< std::endl;
if (parset.get<bool>("io.vtk.write")) {
- std::vector<ScalarVector> stress(assemblers.size());
+ std::vector<ScalarVector> stress(bodyCount);
- for (size_t i=0; i<stress.size(); i++) {
- assemblers[i]->assembleVonMisesStress(body.getYoungModulus(),
- body.getPoissonRatio(),
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = levelContactNetwork.body(i);
+ body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(),
programState.u[i], stress[i]);
}
@@ -479,15 +268,18 @@ int main(int argc, char *argv[]) {
};
report(true);
+
// -------------------
// Set up TNNMG solver
// -------------------
+ auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
+ const auto& totalDirichletNodes = levelContactNetwork.totalDirichletNodes();
using NonlinearFactory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, totalDirichletNodes);
- using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
+ using MyUpdater = Updaters<RateUpdater<Vector, Matrix, typename LevelContactNetwork::Body::Function, dims>,
StateUpdater<ScalarVector, Vector>>;
-
+/*
std::vector<double> L(bodyCount, parset.get<double>("boundary.friction.L"));
std::vector<double> V0(bodyCount, parset.get<double>("boundary.friction.V0"));
MyUpdater current(
@@ -538,8 +330,8 @@ int main(int argc, char *argv[]) {
current.rate_->extractAcceleration(programState.a);
current.state_->extractAlpha(programState.alpha);
- for (size_t i=0; i<deformedGridComplex.size() ; i++) {
- deformedGridComplex.setDeformation(programState.u[i], i);
+ for (size_t i=0; i<bodyCount; i++) {
+ levelContactNetwork.body(i)->setDeformation(programState.u[i]);
}
nBodyAssembler.assembleTransferOperator();
nBodyAssembler.assembleObstacle();
@@ -551,7 +343,7 @@ int main(int argc, char *argv[]) {
break;
}
}
-
+*/
} catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
} catch (std::exception &e) {
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index b850b5cc829734e50342d1e812254778f74ad97b..03c7084ae23e1d8d4a14c76f78504531efe42000 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -11,7 +11,7 @@ vtk.write = false
[problem]
finalTime = 1000 # [s]
-bodyCount = 2
+bodyCount = 3
[body]
bulkModulus = 0.5e5 # [Pa]
diff --git a/src/one-body-problem.cc b/src/one-body-problem.cc
index dedc9162466ba3a15476db11dbbaa9dd7fc120e8..a4548bdf22a8ed436e4815c8762558b614410af0 100644
--- a/src/one-body-problem.cc
+++ b/src/one-body-problem.cc
@@ -46,6 +46,7 @@
#include <dune/fufem/hdf5/file.hh>
#include "assemblers.hh"
+#include "boundarycondition.hh"
#include "diameter.hh"
#include "enumparser.hh"
#include "enums.hh"
diff --git a/src/program_state.hh b/src/program_state.hh
index dfacce35338ca3ab701274a3ff8d2c80619d35d3..ac5c10e322e4223f0423a92dc71e71156ea29f13 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -11,9 +11,10 @@
#include <dune/contact/assemblers/nbodyassembler.hh>
-#include <dune/tectonic/body.hh>
+#include <dune/tectonic/bodydata.hh>
#include "assemblers.hh"
+#include "levelcontactnetwork.hh"
#include "matrices.hh"
#include "spatial-solving/solverfactory.hh"
@@ -44,15 +45,20 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
using ScalarVector = ScalarVectorTEMPLATE;
using BodyState = BodyState<Vector, ScalarVector>;
+private:
+ using LocalVector = typename Vector::block_type;
+ //using LocalMatrix = typename Matrix::block_type;
+ const static int dims = LocalVector::dimension;
+
+public:
ProgramState(const std::vector<size_t>& leafVertexCounts)
: bodyCount_(leafVertexCounts.size()),
- bodies(bodyCount_),
+ bodyStates(bodyCount_),
u(bodyCount_),
v(bodyCount_),
a(bodyCount_),
alpha(bodyCount_),
weightedNormalStress(bodyCount_) {
-
for (size_t i=0; i<bodyCount_; i++) {
size_t leafVertexCount = leafVertexCounts[i];
@@ -62,13 +68,13 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
alpha[i].resize(leafVertexCount),
weightedNormalStress[i].resize(leafVertexCount),
- bodies[i] = new BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
+ bodyStates[i] = new BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
}
}
~ProgramState() {
for (size_t i=0; i<bodyCount_; i++) {
- delete bodies[i];
+ delete bodyStates[i];
}
}
@@ -78,26 +84,14 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Set up initial conditions
- template <class Matrix, class GridView>
- void setupInitialConditions(
- const Dune::ParameterTree& parset,
- const Dune::Contact::NBodyAssembler<typename GridView::Grid, Vector>& nBodyAssembler,
- std::vector<std::function<void(double, Vector &)>> externalForces,
- const Matrices<Matrix>& matrices,
- const std::vector<std::shared_ptr<MyAssembler<GridView, Vector::block_type::dimension>>>& assemblers,
- const std::vector<Dune::BitSetVector<Vector::block_type::dimension>>& dirichletNodes,
- const std::vector<Dune::BitSetVector<Vector::block_type::dimension>>& noNodes,
- const std::vector<BoundaryPatch<GridView>>& frictionalBoundaries,
- const Body<Vector::block_type::dimension>& body) {
-
- using LocalVector = typename Vector::block_type;
- //using LocalMatrix = typename Matrix::block_type;
- auto constexpr dims = LocalVector::dimension;
-
+ 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();
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
- const std::vector<Dune::BitSetVector<dims>>& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+ const Dune::BitSetVector<dims>& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
Dune::ParameterTree const &_localParset) {
@@ -120,6 +114,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// assemble full global contact problem
Matrix bilinearForm;
+
nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
Vector totalRhs;
@@ -128,7 +123,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
Vector totalX;
nBodyAssembler.nodalToTransformed(_x, totalX);
- using LinearFactory = SolverFactory<typename GridView::Grid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+ using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
+ using LinearFactory = SolverFactory<DeformedGridType, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
LinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
auto multigridStep = factory.getStep();
@@ -159,14 +155,14 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
ell0[i].resize(u[i].size());
ell0[i] = 0.0;
- // TODO
- //externalForces[i](relativeTime, ell0[i]);
+
+ levelContactNetwork.body(i)->externalForce(relativeTime, ell0[i]);
}
// Initial displacement: Start from a situation of minimal stress,
// which is automatically attained in the case [v = 0 = a].
// Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
- solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u,
+ solveLinearProblem(levelContactNetwork.totalDirichletNodes(), levelContactNetwork.matrices().elasticity, ell0, u,
parset.sub("u0.solver"));
// Initial acceleration: Computed in agreement with Ma = ell0 - Au
@@ -177,14 +173,25 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
// Initial normal stress
- assemblers[i]->assembleWeightedNormalStress(
- frictionalBoundaries[i], weightedNormalStress[i], body.getYoungModulus(),
- body.getPoissonRatio(), u[i]);
- Dune::MatrixVector::subtractProduct(accelerationRHS[i], *matrices.elasticity[i], u[i]);
+ const auto& body = levelContactNetwork.body(i);
+ std::vector<std::shared_ptr<typename LevelContactNetwork<GridType, dims>::Body::LeafBoundaryCondition>> frictionBoundaryConditions;
+ body->leafBoundaryConditions("friction", frictionBoundaryConditions);
+ for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
+ ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
+
+ body->assembler()->assembleWeightedNormalStress(
+ *frictionBoundaryConditions[j]->boundaryPatch(), frictionBoundaryStress, body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(), u[i]);
+
+ weightedNormalStress[i] += frictionBoundaryStress;
+ }
+
+ Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
}
- solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
+ Dune::BitSetVector<dims> noNodes(levelContactNetwork.totalDirichletNodes().size());
+ solveLinearProblem(noNodes, levelContactNetwork.matrices().mass, accelerationRHS, a,
parset.sub("a0.solver"));
}
@@ -192,7 +199,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
const size_t bodyCount_;
public:
- std::vector<BodyState* > bodies;
+ std::vector<BodyState* > bodyStates;
std::vector<Vector> u;
std::vector<Vector> v;
std::vector<Vector> a;
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 726592df93cd52a2356843051ecee65f5fb33eb8..d5f996fdbb23b0a5073c407d9b2144fee91ba014 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -13,10 +13,11 @@
#include "solverfactory.hh"
-template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
-SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
- const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
- const std::vector<std::vector<Dune::BitSetVector<DeformedGrid::dimension>>>& ignoreNodes)
+
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::SolverFactory(
+ const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
+ const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes)
: nBodyAssembler_(nBodyAssembler),
baseEnergyNorm_(baseSolverStep_),
baseSolver_(&baseSolverStep_,
@@ -34,6 +35,7 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+ /*
// create the transfer operators
const int nCouplings = nBodyAssembler_.nCouplings();
const auto grids = nBodyAssembler_.getGrids();
@@ -41,10 +43,9 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
std::vector<const Dune::BitSetVector<1>*> coarseHasObstacle(nCouplings);
std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
- std::vector<const Matrix*> mortarTransfer(nCouplings);
+ std::vector<const MatrixType*> mortarTransfer(nCouplings);
std::vector<std::array<int,2> > gridIdx(nCouplings);
- /*
for (int j=0; j<nCouplings; j++) {
coarseHasObstacle[j] = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary_[i].getVertices();
fineHasObstacle[j] = nBodyAssembler_.nonmortarBoundary_[j][i+1].getVertices();
@@ -59,8 +60,8 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
nBodyAssembler_.localCoordSystems_[i+1],
coarseHasObstacle, fineHasObstacle,
mortarTransfer,
- gridIdx); */
- }
+ gridIdx);
+ }*/
/*
grids[0], *grids[1], i,
contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
@@ -70,14 +71,14 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
multigridStep_->setTransferOperators(transferOperators_);
}
-template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
-SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::~SolverFactory() {
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::~SolverFactory() {
for (auto &&x : transferOperators_)
delete x;
}
-template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
-auto SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::getStep()
+template <class DeformedGrid, class Nonlinearity, class MatrixType, class VectorType>
+auto SolverFactory<DeformedGrid, Nonlinearity, MatrixType, VectorType>::getStep()
-> std::shared_ptr<Step> {
return multigridStep_;
}
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 932ad8db5256c66ecf2738df778f101e8a1df183..ad4ecbf9b1af5be9334a5af705e077ce31b5c855 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -29,37 +29,36 @@ class SolverFactory {
// const size_t dim = DeformedGrid::dimension;
public:
- using Vector = VectorType;
using Matrix = MatrixType;
-
+ using Vector = VectorType;
using DeformedGrid = DeformedGridTEMPLATE;
using Nonlinearity = NonlinearityTEMPLATE;
public:
- using Step = Dune::Contact::NonSmoothNewtonMGStep<Matrix, Vector>;
+ using Step = Dune::Contact::NonSmoothNewtonMGStep<MatrixType, VectorType>;
- SolverFactory(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
- const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes);
+ SolverFactory(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler,
+ const Dune::BitSetVector<DeformedGrid::dimension>& ignoreNodes);
~SolverFactory();
std::shared_ptr<Step> getStep();
- const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& getNBodyAssembler() const {
+ const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& getNBodyAssembler() const {
return nBodyAssembler_;
}
private:
- const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler_;
+ const Dune::Contact::NBodyAssembler<DeformedGrid, VectorType>& nBodyAssembler_;
- ProjectedBlockGSStep<Matrix, Vector> baseSolverStep_;
- EnergyNorm<Matrix, Vector> baseEnergyNorm_;
- LoopSolver<Vector> baseSolver_;
+ ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep_;
+ EnergyNorm<MatrixType, VectorType> baseEnergyNorm_;
+ LoopSolver<VectorType> baseSolver_;
- ProjectedBlockGSStep<Matrix, Vector> presmoother_;
- ProjectedBlockGSStep<Matrix, Vector> postsmoother_;
+ ProjectedBlockGSStep<MatrixType, VectorType> presmoother_;
+ ProjectedBlockGSStep<MatrixType, VectorType> postsmoother_;
std::shared_ptr<Step> multigridStep_;
- std::vector<Dune::Contact::ContactMGTransfer<Vector>* > transferOperators_;
+ std::vector<Dune::Contact::ContactMGTransfer<VectorType>* > transferOperators_;
};
#endif
diff --git a/src/stackedblocksfactory.cc b/src/stackedblocksfactory.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4c31d3a98c229bea5f6cd4640d1178cd49332b80
--- /dev/null
+++ b/src/stackedblocksfactory.cc
@@ -0,0 +1,189 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "diameter.hh"
+
+#include "multi-body-problem-data/bc.hh"
+#include "multi-body-problem-data/cuboidgeometry.hh"
+#include "multi-body-problem-data/myglobalfrictiondata.hh"
+#include "multi-body-problem-data/weakpatch.hh"
+
+#include "stackedblocksfactory.hh"
+
+template <class GridType, int dims>
+void StackedBlocksFactory<GridType, dims>::setBodies() {
+ // set up cuboid geometries
+
+ double const length = 1.00;
+ double const width = 0.27;
+ double const weakLength = 0.20;
+
+ std::vector<LocalVector> origins(this->bodyCount_);
+ std::vector<LocalVector> lowerWeakOrigins(this->bodyCount_);
+ std::vector<LocalVector> upperWeakOrigins(this->bodyCount_);
+
+#if MY_DIM == 3
+ double const depth = 0.60;
+
+ origins[0] = {0,0,0};
+ lowerWeakOrigins[0] = {0.2,0,0};
+ upperWeakOrigins[0] = {0.2,width,0};
+
+ cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width, depth);
+ for (size_t i=1; i<this->bodyCount_-1; i++) {
+ origins[i] = cuboidGeometries_[i-1]->D;
+ lowerWeakOrigins[i] = {0.6,i*width,0};
+ upperWeakOrigins[i] = {0.6,(i+1)*width,0};
+
+ cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width, depth);
+ }
+ const size_t idx = this->bodyCount_-1;
+ origins[idx] = cuboidGeometries_[idx-1]->D;
+ lowerWeakOrigins[idx] = {0.6,idx*width,0};
+ upperWeakOrigins[idx] = {0.6,(idx+1)*width,0};
+
+ cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width, depth);
+#elif MY_DIM == 2
+ origins[0] = {0,0};
+ lowerWeakOrigins[0] = {0.2,0};
+ upperWeakOrigins[0] = {0.2,width};
+
+ cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width);
+ for (size_t i=1; i<this->bodyCount_-1; i++) {
+ origins[i] = cuboidGeometries_[i-1]->D;
+ lowerWeakOrigins[i] = {0.6,i*width};
+ upperWeakOrigins[i] = {0.6,(i+1)*width};
+
+ cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width);
+ }
+ const size_t idx = this->bodyCount_-1;
+ origins[idx] = cuboidGeometries_[idx-1]->D;
+ lowerWeakOrigins[idx] = {0.6,idx*width};
+ upperWeakOrigins[idx] = {0.6,(idx+1)*width};
+
+ cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width);
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+ // set up reference grids
+ gridConstructor_ = new GridsConstructor<GridType>(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
+ auto lowerWeakPatch = lowerWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
+ auto upperWeakPatch = upperWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
+ getWeakPatch<LocalVector>(this->parset_.sub("boundary.friction.weakening"), cuboidGeometry, *lowerWeakPatch, *upperWeakPatch);
+ refine(*grids[i], *lowerWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
+ refine(*grids[i], *upperWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
+
+ // determine minDiameter and maxDiameter
+ double minDiameter = std::numeric_limits<double>::infinity();
+ double maxDiameter = 0.0;
+ for (auto &&e : elements(grids[i]->leafGridView())) {
+ auto const geometry = e.geometry();
+ auto const diam = diameter(geometry);
+ minDiameter = std::min(minDiameter, diam);
+ maxDiameter = std::max(maxDiameter, diam);
+ }
+ std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
+ std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+ }
+
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ this->bodies_[i] = std::make_shared<typename Base::Body>(bodyData_, grids[i]);
+ }
+}
+
+template <class GridType, int dims>
+void StackedBlocksFactory<GridType, dims>::setCouplings() {
+ const auto deformedGrids = this->levelContactNetwork_.deformedGrids();
+
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& cuboidGeometry = *cuboidGeometries_[i];
+ leafFaces_[i] = std::make_shared<LeafFaces>(this->levelContactNetwork_.leafView(i), cuboidGeometry);
+ levelFaces_[i] = std::make_shared<LevelFaces>(this->levelContactNetwork_.levelView(i), cuboidGeometry);
+ }
+
+ auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+ std::shared_ptr<typename Base::CouplingPair::BackEndType> backend = nullptr;
+ for (size_t i=0; i<this->couplings_.size(); i++) {
+ auto& coupling = this->couplings_[i];
+ coupling = std::make_shared<typename Base::CouplingPair>();
+
+ nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
+ mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+
+ coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::CouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+ }
+}
+
+template <class GridType, int dims>
+void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
+ using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
+ using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+ std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+ 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->levelContactNetwork_.body(i);
+ const auto& leafVertexCount = body->leafVertexCount();
+
+ std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
+
+ // friction boundary
+ if (i<this->bodyCount_-1 && upperWeakPatches_[i]->vertices.size()>0) {
+ std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
+ frictionBoundary->setBoundaryPatch(leafFaces_[i]->upper);
+ frictionBoundary->setWeakeningPatch(upperWeakPatches_[i]);
+ frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], lengthScale));
+ body->addBoundaryCondition(frictionBoundary);
+ }
+ if (i>0 && lowerWeakPatches_[i]->vertices.size()>0) {
+ std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
+ frictionBoundary->setBoundaryPatch(leafFaces_[i]->lower);
+ frictionBoundary->setWeakeningPatch(lowerWeakPatches_[i]);
+ frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *lowerWeakPatches_[i], lengthScale));
+ body->addBoundaryCondition(frictionBoundary);
+ }
+
+ // Neumann boundary
+ std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
+ body->addBoundaryCondition(neumannBoundary);
+
+ // Dirichlet Boundary
+ // identify Dirichlet nodes on leaf level
+ std::shared_ptr<Dune::BitSetVector<dims>> dirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
+ for (int j=0; j<leafVertexCount; j++) {
+ if (leafFaces_[i]->right.containsVertex(j))
+ (*dirichletNodes)[j][0] = true;
+
+ if (leafFaces_[i]->lower.containsVertex(j))
+ (*dirichletNodes)[j][0] = true;
+
+ #if MY_DIM == 3
+ if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
+ dirichletNodes->at(j)[2] = true;
+ #endif
+ }
+
+ std::shared_ptr<LeafBoundaryCondition> dirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ dirichletBoundary->setBoundaryPatch(body->leafView(), dirichletNodes);
+ dirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+ body->addBoundaryCondition(dirichletBoundary);
+ }
+}
+
+#include "stackedblocksfactory_tmpl.cc"
diff --git a/src/stackedblocksfactory.hh b/src/stackedblocksfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..572019d7234ca2c2ecb722a265056e4344f89ed4
--- /dev/null
+++ b/src/stackedblocksfactory.hh
@@ -0,0 +1,77 @@
+#ifndef SRC_STACKEDBLOCKSFACTORY_HH
+#define SRC_STACKEDBLOCKSFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "levelcontactnetworkfactory.hh"
+#include "multi-body-problem-data/mybody.hh"
+#include "multi-body-problem-data/mygrids.hh"
+
+template <class GridType, int dims> class StackedBlocksFactory : public LevelContactNetworkFactory<GridType, dims>{
+private:
+ using Base = LevelContactNetworkFactory<GridType, dims>;
+
+ using LocalVector = typename Base::LevelContactNetwork::LocalVector;
+
+ using DeformedLeafGridView = typename Base::LevelContactNetwork::DeformedLeafGridView;
+ using DeformedLevelGridView = typename Base::LevelContactNetwork::DeformedLevelGridView;
+
+ using LevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
+ using LeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
+
+ using LeafFaces = MyFaces<DeformedLeafGridView>;
+ using LevelFaces = MyFaces<DeformedLevelGridView>;
+
+private:
+ const std::shared_ptr<MyBodyData<dims>> bodyData_; // material properties of bodies
+
+ GridsConstructor<GridType>* 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<ConvexPolyhedron<LocalVector>>> lowerWeakPatches_;
+ std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> upperWeakPatches_;
+
+ std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+ std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+ StackedBlocksFactory(const Dune::ParameterTree& parset) :
+ Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
+ bodyData_(std::make_shared<MyBodyData<dims>>(this->parset_, zenith_())),
+ cuboidGeometries_(this->bodyCount_),
+ leafFaces_(this->bodyCount_),
+ levelFaces_(this->bodyCount_),
+ lowerWeakPatches_(this->bodyCount_),
+ upperWeakPatches_(this->bodyCount_),
+ nonmortarPatch_(this->couplingCount_),
+ mortarPatch_(this->couplingCount_)
+ {}
+
+ ~StackedBlocksFactory() {
+ delete gridConstructor_;
+ }
+
+ void setBodies();
+ void setCouplings();
+ 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/stackedblocksfactory_tmpl.cc b/src/stackedblocksfactory_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..173ea904917e531507cada5f42d8c72d86ea7258
--- /dev/null
+++ b/src/stackedblocksfactory_tmpl.cc
@@ -0,0 +1,7 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "explicitgrid.hh"
+
+template class StackedBlocksFactory<Grid, MY_DIM>;
diff --git a/src/time-stepping/rate.cc b/src/time-stepping/rate.cc
index ce4050c7e0f339564010d0f943902edbb7b9bdb1..b94a2a353cef977101784289a2a88633119a65eb 100644
--- a/src/time-stepping/rate.cc
+++ b/src/time-stepping/rate.cc
@@ -11,7 +11,7 @@ std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
initRateUpdater(Config::scheme scheme,
const std::vector<const Function*>& velocityDirichletFunctions,
const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
- const Matrices<Matrix>& matrices,
+ const Matrices<Matrix,2>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
const std::vector<Vector>& a_initial) {
diff --git a/src/time-stepping/rate.hh b/src/time-stepping/rate.hh
index 5d2ae42809488b41e1021dbc56c29e4313560216..bef6b561b123e012960fe32eed1f72d0fcf0175f 100644
--- a/src/time-stepping/rate.hh
+++ b/src/time-stepping/rate.hh
@@ -10,8 +10,8 @@ template <class Vector, class Matrix, class Function, int dimension>
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
initRateUpdater(Config::scheme scheme,
const std::vector<const Function* >& velocityDirichletFunctions,
- const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
- const Matrices<Matrix>& matrices,
+ const std::vector<std::vector<Dune::BitSetVector<dimension>>>& velocityDirichletNodes,
+ const Matrices<Matrix,2>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
const std::vector<Vector>& a_initial);
diff --git a/src/time-stepping/rate/backward_euler.cc b/src/time-stepping/rate/backward_euler.cc
index 11b7907f7dc46add2283a415c47583ae3bcae7f7..c4845782846180f98dabf1283415bef265ec6ad3 100644
--- a/src/time-stepping/rate/backward_euler.cc
+++ b/src/time-stepping/rate/backward_euler.cc
@@ -5,7 +5,7 @@
template <class Vector, class Matrix, class Function, size_t dim>
BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
- const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
const std::vector<const Function*>& _dirichletFunctions)
diff --git a/src/time-stepping/rate/backward_euler.hh b/src/time-stepping/rate/backward_euler.hh
index 169323d6fcecea70ef0fc9f49efe26053dd44d5e..26e5831d26a184a51d17af729837225294dca122 100644
--- a/src/time-stepping/rate/backward_euler.hh
+++ b/src/time-stepping/rate/backward_euler.hh
@@ -4,7 +4,7 @@
template <class Vector, class Matrix, class Function, size_t dim>
class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
public:
- BackwardEuler(const Matrices<Matrix> &_matrices, const std::vector<Vector> &_u_initial,
+ BackwardEuler(const Matrices<Matrix,2> &_matrices, const std::vector<Vector> &_u_initial,
const std::vector<Vector> &_v_initial, const std::vector<Vector> &_a_initial,
const std::vector<Dune::BitSetVector<dim> > &_dirichletNodes,
const std::vector<const Function*> &_dirichletFunctions);
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index c88136d77457ff1a153daae13b870e54fecee6f8..b3fb46699fb58f959a6e0ef891ffc0cdca3bae06 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -5,7 +5,7 @@
template <class Vector, class Matrix, class Function, size_t dim>
Newmark<Vector, Matrix, Function, dim>::Newmark(
- const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
const std::vector<const Function*>& _dirichletFunctions)
diff --git a/src/time-stepping/rate/newmark.hh b/src/time-stepping/rate/newmark.hh
index 25e014ea2b88821c5945b5f90371ea9435342786..9041e01dc066dde3a03b9285930aaf595598af95 100644
--- a/src/time-stepping/rate/newmark.hh
+++ b/src/time-stepping/rate/newmark.hh
@@ -4,7 +4,7 @@
template <class Vector, class Matrix, class Function, size_t dim>
class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
public:
- Newmark(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ Newmark(const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
const std::vector<const Function*>& _dirichletFunctions);
diff --git a/src/time-stepping/rate/rateupdater.cc b/src/time-stepping/rate/rateupdater.cc
index 6c095235963b39e540e03ccd928c002296b391c8..fb2e22a2d78caa630774df09d9d1aa79cd1b52c2 100644
--- a/src/time-stepping/rate/rateupdater.cc
+++ b/src/time-stepping/rate/rateupdater.cc
@@ -6,7 +6,7 @@
template <class Vector, class Matrix, class Function, size_t dim>
RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
- const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
const std::vector<const Function*>& _dirichletFunctions)
diff --git a/src/time-stepping/rate/rateupdater.hh b/src/time-stepping/rate/rateupdater.hh
index 8b78e1d30cf5b97ebad4db18d454d1e703443651..6a2c73cb1a210c090a257ed136ec6ead3bf9bfba 100644
--- a/src/time-stepping/rate/rateupdater.hh
+++ b/src/time-stepping/rate/rateupdater.hh
@@ -10,7 +10,7 @@
template <class Vector, class Matrix, class Function, size_t dim>
class RateUpdater {
protected:
- RateUpdater(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ RateUpdater(const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
const std::vector<const Function*>& _dirichletFunctions);
@@ -30,7 +30,7 @@ class RateUpdater {
const = 0;
protected:
- const Matrices<Matrix>& matrices;
+ const Matrices<Matrix,2>& matrices;
std::vector<Vector> u, v, a;
const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
const std::vector<const Function*>& dirichletFunctions;
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index c81bc04481353db945e2d92d3a8b283d34851bd5..b71b59eb3e85b3e0255f0bff47bcbd5a92adba37 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -9,7 +9,7 @@ initRateUpdater<Vector, Matrix, Function, MY_DIM>(
Config::scheme scheme,
const std::vector<const Function* >& velocityDirichletFunctions,
const std::vector<Dune::BitSetVector<MY_DIM>>& velocityDirichletNodes,
- const Matrices<Matrix>& matrices,
+ const Matrices<Matrix, 2>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
const std::vector<Vector>& a_initial);