From e495d7ebdf7365cc6fb6b24dd2e410dd53d5ee5b Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 13 Dec 2019 16:45:42 +0100
Subject: [PATCH] major io update
---
dune/tectonic/assemblers.cc | 28 +-
dune/tectonic/assemblers.hh | 7 +-
dune/tectonic/data-structures/body/body.cc | 7 +-
.../data-structures/body/boundarycondition.hh | 21 +-
.../friction/globalfriction.hh | 3 +-
.../friction/globalratestatefriction.hh | 22 +-
.../data-structures/friction/localfriction.hh | 29 +-
.../data-structures/network/contactnetwork.cc | 10 +-
.../data-structures/network/contactnetwork.hh | 1 +
dune/tectonic/factories/CMakeLists.txt | 3 +
.../factories/stackedblocksfactory.cc | 13 +-
dune/tectonic/factories/strikeslipfactory.cc | 275 ++++++++++
dune/tectonic/factories/strikeslipfactory.hh | 87 ++++
.../factories/strikeslipfactory_tmpl.cc | 8 +
dune/tectonic/factories/threeblocksfactory.cc | 4 +-
dune/tectonic/factories/twoblocksfactory.cc | 22 +-
dune/tectonic/io/CMakeLists.txt | 2 +-
dune/tectonic/io/hdf5-bodywriter.hh | 65 ++-
dune/tectonic/io/hdf5-levelwriter.hh | 65 ---
dune/tectonic/io/hdf5-writer.hh | 108 ++++
.../io/hdf5/frictionalboundary-writer.cc | 47 +-
.../io/hdf5/frictionalboundary-writer.hh | 23 +-
.../io/hdf5/frictionalboundary-writer_tmpl.cc | 18 +-
dune/tectonic/io/hdf5/restrict.hh | 20 +-
dune/tectonic/problem-data/CMakeLists.txt | 2 +
dune/tectonic/problem-data/bc.hh | 19 +-
.../tectonic/problem-data/grid/CMakeLists.txt | 6 +
.../problem-data/grid/trianglegeometry.cc | 212 ++++++++
.../problem-data/grid/trianglegeometry.hh | 88 ++++
.../grid/trianglegeometry_tmpl.cc | 8 +
.../problem-data/grid/trianglegrids.cc | 204 ++++++++
.../problem-data/grid/trianglegrids.hh | 77 +++
.../problem-data/grid/trianglegrids_tmpl.cc | 24 +
.../problem-data/segmented-function.hh | 4 +-
dune/tectonic/problem-data/strikeslipbody.hh | 89 ++++
.../contact/dualmortarcoupling.cc | 22 +-
.../spatial-solving/contact/nbodyassembler.cc | 101 ++--
.../spatial-solving/fixedpointiterator.cc | 22 +-
.../spatial-solving/fixedpointiterator.hh | 2 +-
.../spatial-solving/tnnmg/linearization.hh | 4 +-
dune/tectonic/tests/CMakeLists.txt | 3 +
.../tectonic/tests/contacttest/CMakeLists.txt | 33 ++
.../contacttest/staticcontacttest-2D.cfg | 12 +
.../contacttest/staticcontacttest-3D.cfg | 15 +
.../tests/contacttest/staticcontacttest.cc | 313 +++++++++++
.../tests/contacttest/staticcontacttest.cfg | 65 +++
dune/tectonic/tests/hdf5test/CMakeLists.txt | 33 ++
dune/tectonic/tests/hdf5test/hdf5test-2D.cfg | 27 +
dune/tectonic/tests/hdf5test/hdf5test.cc | 202 ++++++++
dune/tectonic/tests/hdf5test/hdf5test.cfg | 105 ++++
.../time-stepping/adaptivetimestepper.cc | 95 +++-
.../state/ageinglawstateupdater.cc | 12 +-
.../state/sliplawstateupdater.cc | 7 +-
dune/tectonic/utils/reductionfactors.hh | 1 +
src/CMakeLists.txt | 1 +
src/foam/foam-2D.cfg | 4 +-
src/foam/foam.cc | 40 +-
src/foam/foam.cfg | 59 ++-
.../multi-body-problem-2D.cfg | 2 +-
src/multi-body-problem/multi-body-problem.cc | 67 ++-
src/multi-body-problem/multi-body-problem.cfg | 12 +-
src/strikeslip/CMakeLists.txt | 43 ++
src/strikeslip/strikeslip-2D.cfg | 27 +
src/strikeslip/strikeslip.cc | 485 ++++++++++++++++++
src/strikeslip/strikeslip.cfg | 112 ++++
65 files changed, 3198 insertions(+), 349 deletions(-)
create mode 100644 dune/tectonic/factories/strikeslipfactory.cc
create mode 100644 dune/tectonic/factories/strikeslipfactory.hh
create mode 100644 dune/tectonic/factories/strikeslipfactory_tmpl.cc
delete mode 100644 dune/tectonic/io/hdf5-levelwriter.hh
create mode 100644 dune/tectonic/io/hdf5-writer.hh
create mode 100644 dune/tectonic/problem-data/grid/trianglegeometry.cc
create mode 100644 dune/tectonic/problem-data/grid/trianglegeometry.hh
create mode 100644 dune/tectonic/problem-data/grid/trianglegeometry_tmpl.cc
create mode 100644 dune/tectonic/problem-data/grid/trianglegrids.cc
create mode 100644 dune/tectonic/problem-data/grid/trianglegrids.hh
create mode 100644 dune/tectonic/problem-data/grid/trianglegrids_tmpl.cc
create mode 100644 dune/tectonic/problem-data/strikeslipbody.hh
create mode 100644 dune/tectonic/tests/contacttest/CMakeLists.txt
create mode 100644 dune/tectonic/tests/contacttest/staticcontacttest-2D.cfg
create mode 100644 dune/tectonic/tests/contacttest/staticcontacttest-3D.cfg
create mode 100644 dune/tectonic/tests/contacttest/staticcontacttest.cc
create mode 100644 dune/tectonic/tests/contacttest/staticcontacttest.cfg
create mode 100644 dune/tectonic/tests/hdf5test/CMakeLists.txt
create mode 100644 dune/tectonic/tests/hdf5test/hdf5test-2D.cfg
create mode 100644 dune/tectonic/tests/hdf5test/hdf5test.cc
create mode 100644 dune/tectonic/tests/hdf5test/hdf5test.cfg
create mode 100644 dune/tectonic/utils/reductionfactors.hh
create mode 100644 src/strikeslip/CMakeLists.txt
create mode 100644 src/strikeslip/strikeslip-2D.cfg
create mode 100644 src/strikeslip/strikeslip.cc
create mode 100644 src/strikeslip/strikeslip.cfg
diff --git a/dune/tectonic/assemblers.cc b/dune/tectonic/assemblers.cc
index 609436ca..8cfc399e 100644
--- a/dune/tectonic/assemblers.cc
+++ b/dune/tectonic/assemblers.cc
@@ -107,13 +107,33 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNeumann(
- BoundaryPatch<GridView> const &neumannBoundary,
- Vector &f,
- Dune::VirtualFunction<double, double> const &neumann,
+ const BoundaryPatch<GridView>& neumannBoundary,
+ const Dune::BitSetVector<dimension>& neumannNodes,
+ Vector& f,
+ const Dune::VirtualFunction<double, double>& neumann,
double relativeTime) const {
+ typename LocalVector::field_type val = 0;
+ neumann.evaluate(relativeTime, val);
+
+ size_t idx = 0;
+ bool neumannIdx = false;
+ for (; idx<neumannNodes.size() && !neumannIdx; idx++) {
+ for (size_t d=0; d<dimension; d++) {
+ if (neumannNodes[idx][d]) {
+ neumannIdx = true;
+ break;
+ }
+ }
+ }
+ idx--;
+
LocalVector localNeumann(0);
- neumann.evaluate(relativeTime, localNeumann[0]);
+ for (size_t i=0; i<localNeumann.size(); i++) {
+ if (neumannNodes[idx][i])
+ localNeumann[i] = val;
+ }
+
NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
localNeumann));
diff --git a/dune/tectonic/assemblers.hh b/dune/tectonic/assemblers.hh
index df0b3cd6..9ec78078 100644
--- a/dune/tectonic/assemblers.hh
+++ b/dune/tectonic/assemblers.hh
@@ -79,9 +79,10 @@ class MyAssembler {
Vector &f) const;
void assembleNeumann(
- BoundaryPatch<GridView> const &neumannBoundary,
- Vector &f,
- Dune::VirtualFunction<double, double> const &neumann,
+ const BoundaryPatch<GridView>& neumannBoundary,
+ const Dune::BitSetVector<dimension>& neumannNodes,
+ Vector& f,
+ const Dune::VirtualFunction<double, double>& neumann,
double relativeTime) const;
void assembleWeightedNormalStress(
diff --git a/dune/tectonic/data-structures/body/body.cc b/dune/tectonic/data-structures/body/body.cc
index 75d92f7d..cd63c1d5 100644
--- a/dune/tectonic/data-structures/body/body.cc
+++ b/dune/tectonic/data-structures/body/body.cc
@@ -34,8 +34,11 @@ LeafBody<GridTEMPLATE, VectorType>::LeafBody(
const auto& leafNeumannCondition = leafNeumannConditions[i];
VectorType b;
- assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(), b, *leafNeumannCondition->boundaryFunction(),
- relativeTime);
+ assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(),
+ *leafNeumannCondition->boundaryNodes(),
+ b,
+ *leafNeumannCondition->boundaryFunction(),
+ relativeTime);
ell += b;
}
};
diff --git a/dune/tectonic/data-structures/body/boundarycondition.hh b/dune/tectonic/data-structures/body/boundarycondition.hh
index 088ebc97..974c5ecb 100644
--- a/dune/tectonic/data-structures/body/boundarycondition.hh
+++ b/dune/tectonic/data-structures/body/boundarycondition.hh
@@ -6,6 +6,7 @@
#include <dune/fufem/boundarypatch.hh>
+#include <dune/tectonic/utils/tobool.hh>
template <class GridView, int dims>
class BoundaryCondition {
@@ -28,9 +29,10 @@ class BoundaryCondition {
BoundaryCondition(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function, const std::string& tag = "") :
tag_(tag),
- boundaryPatch_(patch),
boundaryFunction_(function)
- {}
+ {
+ setBoundaryPatch(patch);
+ }
void setBoundaryPatch(const GridView& gridView, std::shared_ptr<Dune::BitSetVector<dims>> nodes) {
boundaryNodes_ = nodes;
@@ -39,10 +41,21 @@ class BoundaryCondition {
void setBoundaryPatch(std::shared_ptr<BoundaryPatch> patch) {
boundaryPatch_ = patch;
+
+ auto nodes = patch->getVertices();
+ boundaryNodes_ = std::make_shared<Dune::BitSetVector<dims>>(nodes->size(), false);
+ for (size_t i=0; i<nodes->size(); i++) {
+ if (toBool((*nodes)[i])) {
+ for (size_t d=0; d<dims; d++) {
+ (*boundaryNodes_)[i][d] = true;
+ }
+ }
+ }
}
void setBoundaryPatch(const BoundaryPatch& patch) {
- boundaryPatch_ = std::make_shared<BoundaryPatch>(patch);
+ auto patch_ptr = std::make_shared<BoundaryPatch>(patch);
+ setBoundaryPatch(patch_ptr);
}
void setBoundaryFunction(std::shared_ptr<Function> function) {
@@ -50,7 +63,7 @@ class BoundaryCondition {
}
void set(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function) {
- boundaryPatch_ = patch;
+ setBoundaryPatch(patch);
boundaryFunction_ = function;
}
diff --git a/dune/tectonic/data-structures/friction/globalfriction.hh b/dune/tectonic/data-structures/friction/globalfriction.hh
index f3dbabc7..e8796c30 100644
--- a/dune/tectonic/data-structures/friction/globalfriction.hh
+++ b/dune/tectonic/data-structures/friction/globalfriction.hh
@@ -39,8 +39,9 @@ class GlobalFriction {
LocalNonlinearity const virtual &restriction(size_t i) const = 0;
void addHessian(Vector const &v, Matrix &hessian) const {
- for (size_t i = 0; i < v.size(); ++i)
+ for (size_t i = 0; i < v.size(); ++i) {
restriction(i).addHessian(v[i], hessian[i][i]);
+ }
}
void directionalDomain(Vector const &, Vector const &,
diff --git a/dune/tectonic/data-structures/friction/globalratestatefriction.hh b/dune/tectonic/data-structures/friction/globalratestatefriction.hh
index ec42936b..4663d516 100644
--- a/dune/tectonic/data-structures/friction/globalratestatefriction.hh
+++ b/dune/tectonic/data-structures/friction/globalratestatefriction.hh
@@ -38,8 +38,8 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
size_t bodyIndex(const size_t globalIdx) {
size_t i=0;
- for (; i<maxIndex_.size(); i++) {
- if (globalIdx < maxIndex_[i])
+ for (; i<offSet_.size()-1; i++) {
+ if (globalIdx < offSet_[i])
break;
}
return i;
@@ -56,6 +56,10 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
assert(weights.size() == weightedNormalStress.size());
const auto nBodies = weights.size();
+ offSet_.resize(nBodies, 0);
+ for (size_t i=1; i<nBodies; i++) {
+ offSet_[i] = weights[i-1].size();
+ }
std::vector<std::vector<int>> nonmortarBodies(nBodies); // first index body, second index coupling
for (size_t i=0; i<contactCouplings.size(); i++) {
@@ -63,9 +67,6 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
nonmortarBodies[nonmortarIdx].emplace_back(i);
}
- maxIndex_.resize(nBodies);
- size_t globalIdx = 0;
-
for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
const auto& couplingIndices = nonmortarBodies[bodyIdx];
@@ -86,26 +87,23 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
continue;
- localToGlobal_.emplace_back(i);
+ localToGlobal_.emplace_back(offSet_[bodyIdx] + i);
restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
couplings[j]->frictionData()(geoToPoint(it->geometry())));
break;
}
-
- globalIdx++;
}
- maxIndex_[bodyIdx] = globalIdx;
}
}
void updateAlpha(const std::vector<ScalarVector>& alpha) override {
for (size_t j = 0; j < restrictions_.size(); ++j) {
- const auto bodyIdx = bodyIndex(j);
const auto globalDof = localToGlobal_[j];
+ const auto bodyIdx = bodyIndex(globalDof);
size_t bodyDof;
if (bodyIdx>0) {
- bodyDof = globalDof - maxIndex_[bodyIdx-1];
+ bodyDof = globalDof - offSet_[bodyIdx-1];
} else {
bodyDof = globalDof;
}
@@ -125,7 +123,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
private:
std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions_;
- std::vector<size_t> maxIndex_; // max global index per body
+ std::vector<size_t> offSet_; // index off-set by body id
std::vector<size_t> localToGlobal_;
WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction_;
};
diff --git a/dune/tectonic/data-structures/friction/localfriction.hh b/dune/tectonic/data-structures/friction/localfriction.hh
index 67b7f282..461dcea9 100644
--- a/dune/tectonic/data-structures/friction/localfriction.hh
+++ b/dune/tectonic/data-structures/friction/localfriction.hh
@@ -94,7 +94,7 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
//std::cout << A << std::endl;
//std::cout << x << std::endl;
- double const xnorm2 = x.two_norm2();
+ /* double const xnorm2 = x.two_norm2();
double const xnorm = std::sqrt(xnorm2);
if (xnorm2 <= 0.0)
return;
@@ -121,6 +121,33 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
for (size_t k = 0; k < dimension; ++k) {
double const entry = tensorweight * x[k] * x[k];
A[k][k] += entry + idweight;
+ }*/
+
+ const double xnorm = x.two_norm();
+ if (xnorm <= 0.0)
+ return;
+
+ VectorType y = x;
+ y /= xnorm;
+
+ double const H1 = func_.differential(xnorm);
+ double const H2 = func_.second_deriv(xnorm);
+
+ double const tensorweight = (H2 - H1 / xnorm);
+ double const idweight = H1 / xnorm;
+
+ //std::cout << tensorweight << " " << idweight << std::endl;
+
+ for (size_t i = 0; i < dimension; ++i)
+ for (size_t j = 0; j < i; ++j) {
+ double const entry = tensorweight * y[i] * y[j];
+ A[i][j] += entry;
+ A[j][i] += entry;
+ }
+
+ for (size_t k = 0; k < dimension; ++k) {
+ double const entry = tensorweight * y[k] * y[k];
+ A[k][k] += entry + idweight;
}
//std::cout << A << std::endl;
diff --git a/dune/tectonic/data-structures/network/contactnetwork.cc b/dune/tectonic/data-structures/network/contactnetwork.cc
index d9779a74..7e03c549 100644
--- a/dune/tectonic/data-structures/network/contactnetwork.cc
+++ b/dune/tectonic/data-structures/network/contactnetwork.cc
@@ -352,7 +352,7 @@ void ContactNetwork<HostGridType, VectorType>::totalNodes(
const auto& nodes = boundaryConditions[bc]->boundaryNodes();
for (size_t j=0; j<nodes->size(); j++, idx++)
for (int k=0; k<dim; k++)
- totalNodes[idx][k] = (*nodes)[j][k];
+ totalNodes[idx][k] = totalNodes[idx][k] or (*nodes)[j][k];
idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
}
} else {
@@ -421,4 +421,12 @@ void ContactNetwork<HostGridType, VectorType>::frictionNodes(std::vector<const D
}
}
+template <class HostGridType, class VectorType>
+void ContactNetwork<HostGridType, VectorType>::frictionPatches(std::vector<const BoundaryPatch*>& patches) const {
+ patches.resize(nBodies());
+ for (size_t i=0; i<nBodies(); i++) {
+ patches[i] = frictionBoundaries_[i].get();
+ }
+}
+
#include "contactnetwork_tmpl.cc"
diff --git a/dune/tectonic/data-structures/network/contactnetwork.hh b/dune/tectonic/data-structures/network/contactnetwork.hh
index 87b01365..6448fc5f 100644
--- a/dune/tectonic/data-structures/network/contactnetwork.hh
+++ b/dune/tectonic/data-structures/network/contactnetwork.hh
@@ -142,5 +142,6 @@ class ContactNetwork {
const std::string& tag,
BoundaryFunctions& functions) const;
void frictionNodes(std::vector<const Dune::BitSetVector<1>*>& nodes) const;
+ void frictionPatches(std::vector<const BoundaryPatch*>& patches) const;
};
#endif
diff --git a/dune/tectonic/factories/CMakeLists.txt b/dune/tectonic/factories/CMakeLists.txt
index b069d44e..0a757b9d 100644
--- a/dune/tectonic/factories/CMakeLists.txt
+++ b/dune/tectonic/factories/CMakeLists.txt
@@ -5,6 +5,8 @@ add_custom_target(tectonic_dune_factories SOURCES
levelcontactnetworkfactory.hh
stackedblocksfactory.hh
stackedblocksfactory.cc
+ strikeslipfactory.hh
+ strikeslipfactory.cc
threeblocksfactory.hh
threeblocksfactory.cc
twoblocksfactory.hh
@@ -17,6 +19,7 @@ install(FILES
contactnetworkfactory.hh
levelcontactnetworkfactory.hh
stackedblocksfactory.hh
+ strikeslipfactory.hh
threeblocksfactory.hh
twoblocksfactory.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/factories/stackedblocksfactory.cc b/dune/tectonic/factories/stackedblocksfactory.cc
index 6e6f5428..b4665f04 100644
--- a/dune/tectonic/factories/stackedblocksfactory.cc
+++ b/dune/tectonic/factories/stackedblocksfactory.cc
@@ -81,9 +81,8 @@ void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height_i);
-
- cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
+ cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
}
const size_t idx = this->bodyCount_-1;
@@ -145,11 +144,11 @@ void StackedBlocksFactory<HostGridType, VectorType>::setCouplings() {
auto& coupling = this->couplings_[i];
coupling = std::make_shared<typename Base::FrictionCouplingPair>();
- nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
- mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
- weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i]->weakeningRegions()[0]);
+ nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+ mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
+ weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i+1]->weakeningRegions()[0]);
- coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+ coupling->set(i+1, i, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
coupling->setWeakeningPatch(weakPatches_[i]);
coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
}
@@ -183,7 +182,7 @@ void StackedBlocksFactory<HostGridType, VectorType>::setBoundaryConditions() {
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
- std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
const double lengthScale = CuboidGeometry::lengthScale();
diff --git a/dune/tectonic/factories/strikeslipfactory.cc b/dune/tectonic/factories/strikeslipfactory.cc
new file mode 100644
index 00000000..130e6dd2
--- /dev/null
+++ b/dune/tectonic/factories/strikeslipfactory.cc
@@ -0,0 +1,275 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../problem-data/bc.hh"
+#include "../problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "strikeslipfactory.hh"
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBodies() {
+ // set up cuboid geometries
+
+ std::array<double, 2> lengths = {this->parset_.template get<double>("body0.length"), this->parset_.template get<double>("body1.length")};
+ std::array<double, 2> heights = {this->parset_.template get<double>("body0.height"), this->parset_.template get<double>("body1.height")};
+
+ GlobalCoords origin(0);
+
+ const auto& frictionParset = this->parset_.sub("boundary.friction");
+
+#if MY_DIM == 3 // TODO: not implemented
+ std::array<double, 2> depths = {this->parset_.template get<double>("body0.depth"), this->parset_.template get<double>("body1.depth")};
+
+ cuboidGeometries_[0] = std::make_shared<TriangleGeometry>(origins[0], lengths[0], heights[0], depths[0]);
+ cuboidGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
+
+ cuboidGeometries_[1] = std::make_shared<TriangleGeometry>(origins[1], lengths[1], heights[1], depths[1]);
+ cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
+#elif MY_DIM == 2
+ triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origin, lengths[0], heights[0]);
+ triangleGeometries_[0]->addWeakeningPatch(frictionParset, triangleGeometries_[0]->A(), triangleGeometries_[0]->C());
+
+ triangleGeometries_[1] = std::make_shared<TriangleGeometry>(triangleGeometries_[0]->C(), -lengths[1], -heights[1]);
+ triangleGeometries_[1]->addWeakeningPatch(frictionParset, triangleGeometries_[1]->A(), triangleGeometries_[1]->C());
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+ // set up reference grids
+ gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(triangleGeometries_);
+ auto& grids = gridConstructor_->getGrids();
+
+ std::array<double, 2> smallestDiameter = {this->parset_.template get<double>("body0.smallestDiameter"), this->parset_.template get<double>("body1.smallestDiameter")};
+
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& triangleGeometry = *triangleGeometries_[i];
+
+ // define weak patch and refine grid
+ const auto& weakeningRegions = triangleGeometry.weakeningRegions();
+ for (size_t j=0; j<weakeningRegions.size(); j++) {
+ refine(*grids[i], weakeningRegions[j], smallestDiameter[i], TriangleGeometry::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;
+ }
+
+#if MY_DIM == 2
+ bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), 0.0, zenith_());
+ this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
+
+ bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), 0.0, zenith_());
+ this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
+#else
+ bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("gravity"), zenith_());
+ this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
+
+ bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("gravity"), zenith_());
+ this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
+#endif
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
+ const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel()});
+
+ for (size_t l=0; l<=maxLevel; l++) {
+ std::vector<size_t> bodyLevels(2, l);
+ this->contactNetwork_.addLevel(bodyLevels, l);
+ }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& triangleGeometry = *triangleGeometries_[i];
+ leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), triangleGeometry);
+ levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), triangleGeometry);
+ }
+
+ auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+
+ nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->b);
+ mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->b);
+ weakPatches_[0] = std::make_shared<WeakeningRegion>(triangleGeometries_[1]->weakeningRegions()[0]);
+
+ auto& coupling = this->couplings_[0];
+ coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+ coupling->set(1, 0, nonmortarPatch_[0], mortarPatch_[0], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+ coupling->setWeakeningPatch(weakPatches_[0]);
+ coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[0], TriangleGeometry::lengthScale()));
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
+ using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
+
+ using Function = Dune::VirtualFunction<double, double>;
+ std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(-1.0*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
+
+ const double lengthScale = TriangleGeometry::lengthScale();
+
+ // body0
+ const auto& body0 = this->contactNetwork_.body(0);
+ const auto& leafVertexCount0 = body0->nVertices();
+ std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
+ std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
+
+ const auto& gridView0 = body0->gridView();
+ const auto& indexSet0 = gridView0.indexSet();
+
+ for (const auto& e : elements(gridView0)) {
+ for (const auto& is : intersections(gridView0, e)) {
+ if (!is.boundary())
+ continue;
+
+ auto isCenter = is.geometry().center();
+
+ // Dirichlet boundary (c)
+ if (leafFaces_[0]->c.contains(is)) {
+ isCenter -= triangleGeometries_[0]->A();
+
+ if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->length()) > 0.1) {
+ const auto inside = is.inside();
+
+ auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+ int n = refElement.size(is.indexInInside(),1,dim);
+
+ for (int i=0; i<n; i++) {
+ int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+ int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
+ (*zeroDirichletNodes)[globalIdx][1] = true;
+ }
+ }
+ }
+
+ // velocity Dirichlet boundary (a)
+ if (leafFaces_[0]->a.contains(is)) {
+ isCenter -= triangleGeometries_[0]->C();
+
+ if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->height()) > 0.1) {
+ const auto inside = is.inside();
+
+ auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+ int n = refElement.size(is.indexInInside(),1,dim);
+
+ for (int i=0; i<n; i++) {
+ int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+ int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
+ (*velocityDirichletNodes)[globalIdx][0] = true;
+ }
+ }
+ }
+ }
+ }
+
+ std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+
+ std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ zeroDirichletBoundary->setBoundaryPatch(body0->gridView(), zeroDirichletNodes);
+ zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+ body0->addBoundaryCondition(zeroDirichletBoundary);
+
+ std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ velocityDirichletBoundary->setBoundaryPatch(body0->gridView(), velocityDirichletNodes);
+ velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+ body0->addBoundaryCondition(velocityDirichletBoundary);
+
+ // body1
+ const auto& body1 = this->contactNetwork_.body(1);
+ const auto& leafVertexCount1 = body1->nVertices();
+ std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes1 = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+ std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+
+ const auto& gridView1 = body1->gridView();
+ const auto& indexSet1 = gridView1.indexSet();
+
+ for (const auto& e : elements(gridView1)) {
+ for (const auto& is : intersections(gridView1, e)) {
+ if (!is.boundary())
+ continue;
+
+ auto isCenter = is.geometry().center();
+
+ // Neumann boundary, normal load (c)
+ if (leafFaces_[1]->c.contains(is)) {
+ isCenter -= triangleGeometries_[1]->A();
+
+ if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->length()) > 0.1) {
+ const auto inside = is.inside();
+
+ auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+ int n = refElement.size(is.indexInInside(),1,dim);
+
+ for (int i=0; i<n; i++) {
+ int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+ int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
+ (*loadNeumannNodes)[globalIdx][1] = true;
+ }
+ }
+ }
+
+ // Dirichlet boundary (a)
+ if (leafFaces_[1]->a.contains(is)) {
+ isCenter -= triangleGeometries_[1]->C();
+
+ if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->height()) > 0.1) {
+ const auto inside = is.inside();
+
+ auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+ int n = refElement.size(is.indexInInside(),1,dim);
+
+ for (int i=0; i<n; i++) {
+ int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+ int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
+ (*zeroDirichletNodes1)[globalIdx][0] = true;
+ }
+ }
+ }
+ }
+ }
+
+ std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
+
+ auto loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+ loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
+ loadNeumannBoundary->setBoundaryFunction(constantFunction);
+ body1->addBoundaryCondition(loadNeumannBoundary);
+
+ std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary1 = std::make_shared<LeafBoundaryCondition>("dirichlet");
+ zeroDirichletBoundary1->setBoundaryPatch(body1->gridView(), zeroDirichletNodes1);
+ zeroDirichletBoundary1->setBoundaryFunction(zeroFunction);
+ body1->addBoundaryCondition(zeroDirichletBoundary1);
+
+ // body0, body1: natural boundary conditions
+ for (size_t i=0; i<this->bodyCount_; i++) {
+ const auto& body = this->contactNetwork_.body(i);
+
+ // Neumann boundary
+ auto neumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+ auto neumannPatch = std::make_shared<LeafBoundaryPatch>(body->gridView(), true);
+ neumannBoundary->setBoundaryPatch(neumannPatch);
+ neumannBoundary->setBoundaryFunction(neumannFunction);
+ body->addBoundaryCondition(neumannBoundary);
+ }
+}
+
+#include "strikeslipfactory_tmpl.cc"
diff --git a/dune/tectonic/factories/strikeslipfactory.hh b/dune/tectonic/factories/strikeslipfactory.hh
new file mode 100644
index 00000000..0c3d551e
--- /dev/null
+++ b/dune/tectonic/factories/strikeslipfactory.hh
@@ -0,0 +1,87 @@
+#ifndef DUNE_TECTONIC_FACTORIES_STRIKESLIPFACTORY_HH
+#define DUNE_TECTONIC_FACTORIES_STRIKESLIPFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../problem-data/strikeslipbody.hh"
+#include "../problem-data/grid/trianglegrids.hh"
+#include "../problem-data/grid/trianglegeometry.hh"
+
+template <class HostGridType, class VectorType> class StrikeSlipFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+ using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+ using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+ using GlobalCoords = typename ContactNetwork::LocalVector;
+
+ using LeafGridView = typename ContactNetwork::GridView;
+ using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+ using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+ using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+ using LeafFaces = MyFaces<LeafGridView>;
+ using LevelFaces = MyFaces<LevelGridView>;
+
+ using TriangleGeometry= TriangleGeometry<typename GlobalCoords::field_type>;
+ using WeakeningRegion = typename TriangleGeometry::WeakeningRegion;
+
+ static const int dim = ContactNetwork::dim;
+
+ std::vector<std::shared_ptr<MyBodyData<dim>>> bodyData_; // material properties of bodies
+
+ std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
+
+ std::vector<std::shared_ptr<TriangleGeometry>> triangleGeometries_;
+
+ std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+ std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+ std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+ std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+ std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+ StrikeSlipFactory(const Dune::ParameterTree& parset) :
+ Base(parset, 2, 1),
+ bodyData_(2),
+ triangleGeometries_(2),
+ leafFaces_(2),
+ levelFaces_(2),
+ weakPatches_(2),
+ nonmortarPatch_(1),
+ mortarPatch_(1)
+ {}
+
+ void setBodies();
+ void setLevelBodies();
+ void setCouplings();
+ void setLevelCouplings() {}
+ void setBoundaryConditions();
+
+ auto& weakPatches() {
+ return weakPatches_;
+ }
+
+private:
+ static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+ #if MY_DIM == 2
+ return {0, 0};
+ #else
+ return {0, 0, 1};
+ #endif
+ }
+};
+#endif
+
+
diff --git a/dune/tectonic/factories/strikeslipfactory_tmpl.cc b/dune/tectonic/factories/strikeslipfactory_tmpl.cc
new file mode 100644
index 00000000..ee3446fb
--- /dev/null
+++ b/dune/tectonic/factories/strikeslipfactory_tmpl.cc
@@ -0,0 +1,8 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class StrikeSlipFactory<Grid, Vector>;
diff --git a/dune/tectonic/factories/threeblocksfactory.cc b/dune/tectonic/factories/threeblocksfactory.cc
index 6c395565..2f6fd5cf 100644
--- a/dune/tectonic/factories/threeblocksfactory.cc
+++ b/dune/tectonic/factories/threeblocksfactory.cc
@@ -17,7 +17,7 @@ template <class HostGridType, class VectorTEMPLATE>
void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
// set up cuboid geometries
- std::array<double, 3> lengths = {3.0, 1.0, 1.0};
+ std::array<double, 3> lengths = {2.0, 1.0, 1.0};
std::array<double, 3> heights = {1.0, 1.0, 1.0};
std::array<GlobalCoords, 3> origins;
@@ -147,7 +147,7 @@ void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
- std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
const double lengthScale = CuboidGeometry::lengthScale();
diff --git a/dune/tectonic/factories/twoblocksfactory.cc b/dune/tectonic/factories/twoblocksfactory.cc
index 39c47cd7..44070552 100644
--- a/dune/tectonic/factories/twoblocksfactory.cc
+++ b/dune/tectonic/factories/twoblocksfactory.cc
@@ -121,7 +121,7 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
- std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
const double lengthScale = CuboidGeometry::lengthScale();
@@ -148,7 +148,7 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
const auto& body1 = this->contactNetwork_.body(1);
const auto& leafVertexCount1 = body1->nVertices();
std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
- for (int j=0; j<leafVertexCount1; j++) {
+ for (size_t j=0; j<leafVertexCount1; j++) {
if (leafFaces_[1]->upper.containsVertex(j))
(*velocityDirichletNodes)[j][0] = true;
@@ -167,17 +167,21 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
// body0, body1: natural boundary conditions
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& body = this->contactNetwork_.body(i);
- const auto& leafVertexCount = body->nVertices();
// Neumann boundary
- std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
+ auto neumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+ auto neumannPatch = std::make_shared<LeafBoundaryPatch>(body->gridView(), true);
+ neumannBoundary->setBoundaryPatch(neumannPatch);
+ neumannBoundary->setBoundaryFunction(neumannFunction);
body->addBoundaryCondition(neumannBoundary);
}
// body1: Neumann boundary (upper)
// normal load
- /*std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
- for (int j=0; j<leafVertexCount1; j++) {
+ std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
+
+ std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+ for (size_t j=0; j<leafVertexCount1; j++) {
if (leafFaces_[1]->upper.containsVertex(j))
(*loadNeumannNodes)[j][1] = true;
@@ -187,12 +191,10 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
#endif
}
- std::shared_ptr<LeafBoundaryCondition> loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
- std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
-
+ auto loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
loadNeumannBoundary->setBoundaryFunction(constantFunction);
- body1->addBoundaryCondition(loadNeumannBoundary);*/
+ body1->addBoundaryCondition(loadNeumannBoundary);
}
#include "twoblocksfactory_tmpl.cc"
diff --git a/dune/tectonic/io/CMakeLists.txt b/dune/tectonic/io/CMakeLists.txt
index c1aada25..07045307 100644
--- a/dune/tectonic/io/CMakeLists.txt
+++ b/dune/tectonic/io/CMakeLists.txt
@@ -2,7 +2,7 @@ add_subdirectory("hdf5")
add_custom_target(tectonic_dune_io SOURCES
hdf5-bodywriter.hh
- hdf5-levelwriter.hh
+ hdf5-writer.hh
uniform-grid-writer.cc
vtk.hh
vtk.cc
diff --git a/dune/tectonic/io/hdf5-bodywriter.hh b/dune/tectonic/io/hdf5-bodywriter.hh
index dfb2b3aa..91dd5ede 100644
--- a/dune/tectonic/io/hdf5-bodywriter.hh
+++ b/dune/tectonic/io/hdf5-bodywriter.hh
@@ -6,16 +6,18 @@
#include <dune/fufem/hdf5/file.hh>
#include "hdf5/frictionalboundary-writer.hh"
-#include "hdf5/patchinfo-writer.hh"
+//#include "hdf5/patchinfo-writer.hh"
template <class ProgramState, class VertexBasis, class GridView>
class HDF5BodyWriter {
private:
- using PatchInfoWriter = PatchInfoWriter<ProgramState, VertexBasis, GridView>;
- using FrictionalBoundaryWriter = FrictionalBoundaryWriter<ProgramState, GridView>;
+ using Vector = typename ProgramState::Vector;
+ using ScalarVector = typename ProgramState::ScalarVector;
+ //using PatchInfoWriter = PatchInfoWriter<ProgramState, VertexBasis, GridView>;
+
public:
using VertexCoordinates = typename ProgramState::Vector;
- using FrictionPatches = std::vector<const BoundaryPatch<GridView>* >;
+ using Patch = typename FrictionalBoundaryWriter<GridView>::Patch;
using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
using LocalVector = typename VertexCoordinates::block_type;
@@ -23,50 +25,45 @@ class HDF5BodyWriter {
//friend class HDF5LevelWriter<ProgramState, VertexBasis, GridView>;
HDF5BodyWriter(
- HDF5::Grouplike &file,
+ HDF5::File& file,
+ const size_t bodyID,
const VertexCoordinates& vertexCoordinates,
const VertexBasis& vertexBasis,
- const FrictionPatches& frictionPatches,
+ const Patch& frictionPatch,
const WeakPatches& weakPatches) :
- patchCount_(frictionPatches.size()),
- file_(file),
-#if MY_DIM == 3
- patchInfoWriters_(patchCount_),
-#endif
- frictionBoundaryWriters_(patchCount_) {
+ id_(bodyID),
+ /*#if MY_DIM == 3
+ patchInfoWriters_(patchCount_),
+ #endif*/
+ group_(file, "body"+std::to_string(id_)) {
- assert(patchCount_ == weakPatches.size());
+ frictionBoundaryWriter_ = std::make_unique<FrictionalBoundaryWriter<GridView>>(group_, vertexCoordinates, frictionPatch);
- for (size_t i=0; i<patchCount_; i++) {
-#if MY_DIM == 3
+ /*#if MY_DIM == 3
patchInfoWriters_[i] = std::make_unique<PatchInfoWriter>(file_, vertexBasis, *frictionPatches[i], *weakPatches[i], i);
-#endif
- frictionBoundaryWriters_[i] = std::make_unique<FrictionalBoundaryWriter>(file_, vertexCoordinates, *frictionPatches[i], i);
- }
+ #endif*/
}
- template <class Friction>
- void reportSolution(ProgramState const &programState,
- // for the friction coefficient
- Friction& friction) {
+ void reportSolution(const ProgramState& programState, const Vector& v, const ScalarVector& frictionCoeff) {
+ frictionBoundaryWriter_->write(programState.timeStep, programState.u[id_], v, programState.alpha[id_], frictionCoeff);
- for (size_t i=0; i<patchCount_; i++) {
-#if MY_DIM == 3
- patchInfoWriters_[i]->write(programState);
-#endif
- frictionBoundaryWriters_[i]->write(programState, friction);
- }
+ /*#if MY_DIM == 3
+ patchInfoWriters_[i]->write(programState);
+ #endif*/
}
-private:
- const size_t patchCount_;
+ auto id() const {
+ return id_;
+ }
- HDF5::Grouplike &file_;
+private:
+ const size_t id_;
+ HDF5::Group group_;
-#if MY_DIM == 3
+/*#if MY_DIM == 3
std::vector<std::unique_ptr<PatchInfoWriter>> patchInfoWriters_;
-#endif
+#endif*/
- std::vector<std::unique_ptr<FrictionalBoundaryWriter>> frictionBoundaryWriters_;
+ std::unique_ptr<FrictionalBoundaryWriter<GridView>> frictionBoundaryWriter_;
};
#endif
diff --git a/dune/tectonic/io/hdf5-levelwriter.hh b/dune/tectonic/io/hdf5-levelwriter.hh
deleted file mode 100644
index 4edabd8b..00000000
--- a/dune/tectonic/io/hdf5-levelwriter.hh
+++ /dev/null
@@ -1,65 +0,0 @@
-#ifndef SRC_HDF5_LEVELWRITER_HH
-#define SRC_HDF5_LEVELWRITER_HH
-
-#include <dune/fufem/hdf5/file.hh>
-
-#include "hdf5/iteration-writer.hh"
-#include "hdf5/time-writer.hh"
-#include "hdf5-bodywriter.hh"
-
-template <class ProgramState, class VertexBasis, class GridView>
-class HDF5LevelWriter {
-public:
- using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
- using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
- using VertexBases = std::vector<const VertexBasis*>;
- using FrictionPatches = std::vector<typename HDF5BodyWriter::FrictionPatches>;
- using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
-
- //friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
-
- HDF5LevelWriter(HDF5::Grouplike &file,
- const VertexCoordinates& vertexCoordinates,
- const VertexBases& vertexBases,
- const FrictionPatches& frictionBoundaries,
- const WeakPatches& weakPatches)
- : bodyCount_(vertexCoordinates.size()),
- file_(file),
- iterationWriter_(file_),
- timeWriter_(file_),
- bodyWriters_(bodyCount_) {
-
- for (size_t i=0; i<bodyCount_; i++) {
- bodyWriters_[i] = std::make_unique<HDF5BodyWriter>(file_, vertexCoordinates[i], *vertexBases[i], frictionBoundaries[i], weakPatches[i]);
- }
- }
-
- template <class Friction>
- void reportSolution(
- const ProgramState& programState,
- // for the friction coefficient
- Friction& friction) {
- timeWriter_.write(programState);
-
- for (size_t i=0; i<bodyCount_; i++) {
- bodyWriters_[i]->reportSolution(programState, friction); //TODO
- }
- }
-
- void reportIterations(
- const ProgramState& programState,
- const IterationRegister& iterationCount) {
- iterationWriter_.write(programState.timeStep, iterationCount);
- }
-
-private:
- const size_t bodyCount_;
-
- HDF5::Grouplike &file_;
-
- IterationWriter iterationWriter_;
- TimeWriter<ProgramState> timeWriter_;
-
- std::vector<std::unique_ptr<HDF5BodyWriter>> bodyWriters_;
-};
-#endif
diff --git a/dune/tectonic/io/hdf5-writer.hh b/dune/tectonic/io/hdf5-writer.hh
new file mode 100644
index 00000000..2a8ecc2f
--- /dev/null
+++ b/dune/tectonic/io/hdf5-writer.hh
@@ -0,0 +1,108 @@
+#ifndef SRC_HDF5_WRITER_HH
+#define SRC_HDF5_WRITER_HH
+
+#include <dune/fufem/hdf5/file.hh>
+
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/time-writer.hh"
+#include "hdf5-bodywriter.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5Writer {
+public:
+ using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
+ using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
+ using VertexBases = std::vector<const VertexBasis*>;
+ using FrictionPatches = std::vector<const typename HDF5BodyWriter::Patch*>;
+ using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
+
+ //friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
+
+ HDF5Writer(HDF5::File& file,
+ const VertexCoordinates& vertexCoordinates,
+ const VertexBases& vertexBases,
+ const FrictionPatches& frictionPatches,
+ const WeakPatches& weakPatches)
+ : file_(file),
+ iterationWriter_(file_),
+ timeWriter_(file_),
+ frictionPatches_(frictionPatches) {
+
+ for (size_t i=0; i<frictionPatches_.size(); i++) {
+ if (frictionPatches_[i]->numVertices() > 0)
+ bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i], weakPatches[i]));
+ }
+ }
+
+ template <class ContactNetwork, class Friction>
+ void reportSolution(const ProgramState& programState, const ContactNetwork& contactNetwork, Friction& friction) {
+
+ timeWriter_.write(programState);
+
+ friction.updateAlpha(programState.alpha);
+
+ // extract relative velocities
+ using Vector = typename ProgramState::Vector;
+ Vector mortarV;
+ contactNetwork.nBodyAssembler().nodalToTransformed(programState.v, mortarV);
+
+ std::vector<Vector> v_rel;
+ split(mortarV, v_rel);
+
+ using ScalarVector = typename ProgramState::ScalarVector;
+ const auto frictionCoeff = friction.coefficientOfFriction(mortarV);
+
+ double norm = 0;
+ const auto& bodyNodes = *frictionPatches_[0]->getVertices();
+ for (size_t i=bodyNodes.size(); i<frictionCoeff.size(); i++) {
+ norm += frictionCoeff[i].two_norm();
+ }
+ std::cout << std::setprecision(10) << "friction coefficients norm: " << norm << std::endl;
+
+ std::vector<ScalarVector> splitCoeff;
+ split(frictionCoeff, splitCoeff);
+
+ //print(v_rel, "v_rel: ");
+
+ for (size_t i=0; i<bodyWriters_.size(); i++) {
+ auto bodyID = bodyWriters_[i]->id();
+ bodyWriters_[i]->reportSolution(programState, v_rel[bodyID], splitCoeff[bodyID]);
+ }
+ }
+
+ void reportIterations(const ProgramState& programState, const IterationRegister& iterationCount) {
+ iterationWriter_.write(programState.timeStep, iterationCount);
+ }
+
+ template <class VectorType>
+ void split(const VectorType& v, std::vector<VectorType>& splitV) const {
+
+ const size_t nBodies = frictionPatches_.size();
+
+ size_t globalIdx = 0;
+ splitV.resize(nBodies);
+ for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+ const auto& bodyNodes = *frictionPatches_[bodyIdx]->getVertices();
+
+ auto& splitV_ = splitV[bodyIdx];
+ splitV_.resize(bodyNodes.size());
+
+ for (size_t i=0; i<splitV_.size(); i++) {
+ if (toBool(bodyNodes[i])) {
+ splitV_[i] = v[globalIdx];
+ }
+ globalIdx++;
+ }
+ }
+ }
+
+private:
+ HDF5::File& file_;
+
+ IterationWriter iterationWriter_;
+ TimeWriter<ProgramState> timeWriter_;
+
+ const FrictionPatches& frictionPatches_;
+ std::vector<std::unique_ptr<HDF5BodyWriter>> bodyWriters_;
+};
+#endif
diff --git a/dune/tectonic/io/hdf5/frictionalboundary-writer.cc b/dune/tectonic/io/hdf5/frictionalboundary-writer.cc
index 0cb5e01a..6bcfa353 100644
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer.cc
+++ b/dune/tectonic/io/hdf5/frictionalboundary-writer.cc
@@ -7,12 +7,12 @@
#include "frictionalboundary-writer.hh"
#include "restrict.hh"
-template <class ProgramState, class GridView>
-FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
- HDF5::Grouplike &file, Vector const &vertexCoordinates,
- Patch const &frictionalBoundary, size_t const id)
- : id_(id),
- group_(file, "frictionalBoundary"+std::to_string(id_)),
+template <class GridView>
+template <class Vector>
+FrictionalBoundaryWriter<GridView>::FrictionalBoundaryWriter(
+ HDF5::Group& group, const Vector& vertexCoordinates,
+ const Patch& frictionalBoundary)
+ : group_(group),
frictionalBoundary_(frictionalBoundary),
frictionalBoundaryDisplacementWriter_(group_, "displacement",
frictionalBoundary.numVertices(),
@@ -32,33 +32,22 @@ FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
}
-template <class ProgramState, class GridView>
-template <class Friction>
-void FrictionalBoundaryWriter<ProgramState, GridView>::write(
- ProgramState const &programState, Friction &friction) {
+template <class GridView>
+template <class Vector, class ScalarVector>
+void FrictionalBoundaryWriter<GridView>::write(
+ const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff) {
- auto const frictionalBoundaryDisplacements =
- restrictToSurface(programState.u[id_], frictionalBoundary_);
- addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
- frictionalBoundaryDisplacements);
+ auto const frictionalBoundaryDisplacements = restrictToSurface(u, frictionalBoundary_);
+ addEntry(frictionalBoundaryDisplacementWriter_, timeStep, frictionalBoundaryDisplacements);
- auto const frictionalBoundaryVelocities =
- restrictToSurface(programState.v[id_], frictionalBoundary_);
- addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
- frictionalBoundaryVelocities);
+ auto const frictionalBoundaryVelocities = restrictToSurface(v, frictionalBoundary_);
+ addEntry(frictionalBoundaryVelocityWriter_, timeStep, frictionalBoundaryVelocities);
- auto const frictionalBoundaryStates =
- restrictToSurface(programState.alpha[id_], frictionalBoundary_);
- addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
- frictionalBoundaryStates);
+ auto const frictionalBoundaryStates = restrictToSurface(alpha, frictionalBoundary_);
+ addEntry(frictionalBoundaryStateWriter_, timeStep, frictionalBoundaryStates);
- //TODO: needs to transform programState.v to relative velocities with nBodyAssembler
- /*friction.updateAlpha(programState.alpha);
- auto const c = friction.coefficientOfFriction(programState.v[id_]);
- auto const frictionalBoundaryCoefficient =
- restrictToSurface(c, frictionalBoundary_);
- addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
- frictionalBoundaryCoefficient);*/
+ auto const frictionalBoundaryCoefficient = restrictToSurface(frictionCoeff, frictionalBoundary_);
+ addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
}
#include "frictionalboundary-writer_tmpl.cc"
diff --git a/dune/tectonic/io/hdf5/frictionalboundary-writer.hh b/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
index 09774013..a7ca0855 100644
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
+++ b/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
@@ -1,28 +1,25 @@
#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
-#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
-template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
- using ScalarVector = typename ProgramState::ScalarVector;
- using Vector = typename ProgramState::Vector;
+template <class GridView> class FrictionalBoundaryWriter {
+public:
using Patch = BoundaryPatch<GridView>;
-public:
- FrictionalBoundaryWriter(HDF5::Grouplike& file, Vector const &vertexCoordinates,
- Patch const &frictionalBoundary, size_t const id);
+ template <class Vector>
+ FrictionalBoundaryWriter(HDF5::Group& group, const Vector& vertexCoordinates,
+ const Patch& frictionalBoundary);
- template <class Friction>
- void write(ProgramState const &programState, Friction &friction);
+ template <class Vector, class ScalarVector>
+ void write(const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff);
private:
- size_t const id_;
-
- HDF5::Group group_;
+ HDF5::Group& group_;
- Patch const &frictionalBoundary_;
+ const Patch& frictionalBoundary_;
HDF5::SequenceIO<2> frictionalBoundaryDisplacementWriter_;
HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
diff --git a/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc b/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc
index 883c59ba..51633c69 100644
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc
+++ b/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc
@@ -1,13 +1,19 @@
#include "../../explicitvectors.hh"
#include "../../explicitgrid.hh"
-#include "../../data-structures/program_state.hh"
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using MyFriction = GlobalFriction<Matrix, Vector>;
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
-template class FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>;
-template void FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>::write(
- MyProgramState const &programState, MyFriction &friction);
+//#include "../../data-structures/program_state.hh"
+
+//using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class FrictionalBoundaryWriter<DefLeafGridView>;
+
+template FrictionalBoundaryWriter<DefLeafGridView>::FrictionalBoundaryWriter(
+ HDF5::Group& group, const Vector& vertexCoordinates, const BoundaryPatch<DefLeafGridView>& frictionalBoundary);
+template void FrictionalBoundaryWriter<DefLeafGridView>::write(
+ const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff);
diff --git a/dune/tectonic/io/hdf5/restrict.hh b/dune/tectonic/io/hdf5/restrict.hh
index 78639f1b..efcdaffc 100644
--- a/dune/tectonic/io/hdf5/restrict.hh
+++ b/dune/tectonic/io/hdf5/restrict.hh
@@ -8,16 +8,16 @@
#include "../../utils/tobool.hh"
template <class Vector, class Patch>
-Vector restrictToSurface(Vector const &v1, Patch const &patch) {
- auto const &vertices = *patch.getVertices();
- assert(vertices.size() == v1.size());
+Vector restrictToSurface(const Vector& v1, const Patch& patch) {
+ auto const &vertices = *patch.getVertices();
+ assert(vertices.size() == v1.size());
- Vector ret(vertices.count());
- auto target = ret.begin();
- for (size_t i = 0; i < v1.size(); ++i)
- if (toBool(vertices[i]))
- *(target++) = v1[i];
- assert(target == ret.end());
- return ret;
+ Vector ret(vertices.count());
+ auto target = ret.begin();
+ for (size_t i = 0; i < v1.size(); ++i)
+ if (toBool(vertices[i]))
+ *(target++) = v1[i];
+ assert(target == ret.end());
+ return ret;
}
#endif
diff --git a/dune/tectonic/problem-data/CMakeLists.txt b/dune/tectonic/problem-data/CMakeLists.txt
index 8b975765..e59f2e04 100644
--- a/dune/tectonic/problem-data/CMakeLists.txt
+++ b/dune/tectonic/problem-data/CMakeLists.txt
@@ -8,6 +8,7 @@ add_custom_target(tectonic_dune_problem-data SOURCES
myglobalfrictiondata.hh
patchfunction.hh
segmented-function.hh
+ strikeslipbody.hh
)
#install headers
@@ -19,4 +20,5 @@ install(FILES
myglobalfrictiondata.hh
patchfunction.hh
segmented-function.hh
+ strikeslipbody.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/problem-data/bc.hh b/dune/tectonic/problem-data/bc.hh
index 2d478035..ac82e47f 100644
--- a/dune/tectonic/problem-data/bc.hh
+++ b/dune/tectonic/problem-data/bc.hh
@@ -5,9 +5,14 @@
class VelocityDirichletCondition
: public Dune::VirtualFunction<double, double> {
+
+public:
+ VelocityDirichletCondition(const double finalVelocity = 5e-5) :
+ finalVelocity_(finalVelocity)
+ {}
+
void evaluate(double const &relativeTime, double &y) const {
// Assumed to vanish at time zero
- double const finalVelocity = 5e-5;
//std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
@@ -17,15 +22,21 @@ class VelocityDirichletCondition
std::cout << "- final velocity reached" << std::endl;*/
y = (relativeTime <= 0.1)
- ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
- : finalVelocity;
+ ? finalVelocity_ * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
+ : finalVelocity_;
+
+ //y = relativeTime * finalVelocity;
+
}
+
+private:
+ const double finalVelocity_;
};
class NeumannCondition : public Dune::VirtualFunction<double, double> {
public:
NeumannCondition(double c = 0.0) : c_(c) {}
- void evaluate(double const &relativeTime, double &y) const { y = c_; }
+ void evaluate(double const &relativeTime, double &y) const { y = -c_; }
private:
double c_;
diff --git a/dune/tectonic/problem-data/grid/CMakeLists.txt b/dune/tectonic/problem-data/grid/CMakeLists.txt
index e0fd59ba..b7abfed5 100644
--- a/dune/tectonic/problem-data/grid/CMakeLists.txt
+++ b/dune/tectonic/problem-data/grid/CMakeLists.txt
@@ -11,6 +11,10 @@ add_custom_target(tectonic_dune_problem-data_grid SOURCES
mygrids.cc
simplexmanager.hh
simplexmanager.cc
+ trianglegeometry.hh
+ trianglegeometry.cc
+ trianglegrids.hh
+ trianglegrids.cc
)
#install headers
@@ -22,4 +26,6 @@ install(FILES
gridconstructor.hh
mygrids.hh
simplexmanager.hh
+ trianglegeometry.hh
+ trianglegrids.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/problem-data/grid/trianglegeometry.cc b/dune/tectonic/problem-data/grid/trianglegeometry.cc
new file mode 100644
index 00000000..7c24275f
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/trianglegeometry.cc
@@ -0,0 +1,212 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+#include "trianglegeometry.hh"
+
+#if MY_DIM == 3
+template <class ctype>
+TriangleGeometry<ctype>::TriangleGeometry(const GlobalCoords& origin,
+ const double length, const double height, const double depth) :
+ length_(length*lengthScale()),
+ height_(height*lengthScale()),
+ depth_(depth*lengthScale()),
+ A_(origin),
+ B_({origin[0]+length_, origin[1], 0}),
+ C_({origin[0]+length_, origin[1]+height_, 0})
+ {}
+#else
+
+template <class ctype>
+TriangleGeometry<ctype>::TriangleGeometry(const GlobalCoords& origin,
+ const double length, const double height) :
+ length_(length*lengthScale()),
+ height_(height*lengthScale()),
+ A_(origin),
+ B_({origin[0]+length_, origin[1]}),
+ C_({origin[0]+length_, origin[1]+height_})
+ {}
+#endif
+
+template <class ctype>
+void TriangleGeometry<ctype>::addWeakeningRegion(const WeakeningRegion& weakeningRegion) {
+ weakeningRegions_.emplace_back(weakeningRegion);
+}
+
+template <class ctype>
+void TriangleGeometry<ctype>::addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1) {
+ WeakeningRegion weakPatch;
+
+#if MY_DIM == 3 // TODO: Does not work yet
+ if (vertex0 != vertex1) {
+ weakPatch.vertices.resize(4);
+ weakPatch.vertices[0] = weakPatch.vertices[2] = vertex0;
+ weakPatch.vertices[1] = weakPatch.vertices[3] = vertex1;
+
+ for (size_t k = 0; k < 2; ++k) {
+ weakPatch.vertices[k][2] = -depth_ / 2.0;
+ weakPatch.vertices[k + 2][2] = depth_ / 2.0;
+ }
+ switch (parset.get<Config::PatchType>("patchType")) {
+ case Config::Rectangular:
+ break;
+ case Config::Trapezoidal:
+ weakPatch.vertices[1][0] += 0.05 * lengthScale();
+ weakPatch.vertices[3][0] -= 0.05 * lengthScale();
+ break;
+ default:
+ assert(false);
+ }
+ addWeakeningRegion(weakPatch);
+ }
+#else
+ if (vertex0 != vertex1) {
+ weakPatch.vertices.resize(2);
+ weakPatch.vertices[0] = vertex0;
+ weakPatch.vertices[1] = vertex1;
+ addWeakeningRegion(weakPatch);
+ }
+#endif
+}
+
+template <class ctype>
+void TriangleGeometry<ctype>::write() const {
+ std::fstream writer("geometry", std::fstream::out);
+ writer << "A = " << A_ << std::endl;
+ writer << "B = " << B_ << std::endl;
+ writer << "C = " << C_ << std::endl;
+}
+
+/*
+template <class ctype>
+void CuboidGeometry<ctype>::render() const {
+#ifdef HAVE_CAIROMM
+ std::string const filename = "geometry.png";
+ double const width = 600;
+ double const height = 400;
+ double const widthScale = 400;
+ double const heightScale = 400;
+
+ auto surface =
+ Cairo::ImageSurface::create(Cairo::FORMAT_ARGB32, width, height);
+ auto cr = Cairo::Context::create(surface);
+
+ auto const setRGBColor = [&](int colour) {
+ cr->set_source_rgb(((colour & 0xFF0000) >> 16) / 255.0,
+ ((colour & 0x00FF00) >> 8) / 255.0,
+ ((colour & 0x0000FF) >> 0) / 255.0);
+ };
+ auto const moveTo = [&](LocalVector2D const &v) { cr->move_to(v[0], -v[1]); };
+ auto const lineTo = [&](LocalVector2D const &v) { cr->line_to(v[0], -v[1]); };
+
+ cr->scale(widthScale, heightScale);
+ cr->translate(0.1, 0.1);
+ cr->set_line_width(0.0025);
+
+ // triangle
+ {
+ moveTo(reference::A);
+ lineTo(reference::B);
+ lineTo(reference::C);
+ cr->close_path();
+ cr->stroke();
+ }
+
+ // dashed lines
+ {
+ cr->save();
+ std::vector<double> dashPattern = { 0.005 };
+ cr->set_dash(dashPattern, 0);
+ moveTo(reference::Z);
+ lineTo(reference::Y);
+ moveTo(reference::U);
+ lineTo(reference::X);
+ cr->stroke();
+ cr->restore();
+ }
+
+ // fill viscoelastic region
+ {
+ cr->save();
+ setRGBColor(0x0097E0);
+ moveTo(reference::B);
+ lineTo(reference::K);
+ lineTo(reference::M);
+ cr->fill();
+ cr->restore();
+ }
+
+ // mark weakening region
+ {
+ cr->save();
+ setRGBColor(0x7AD3FF);
+ cr->set_line_width(0.005);
+ moveTo(reference::X);
+ lineTo(reference::Y);
+ cr->stroke();
+ cr->restore();
+ }
+
+ // mark points
+ {
+ auto const drawCircle = [&](LocalVector2D const &v) {
+ cr->arc(v[0], -v[1], 0.0075, -M_PI, M_PI); // x,y,radius,angle1,angle2
+ cr->fill();
+ };
+
+ cr->save();
+ setRGBColor(0x002F47);
+ drawCircle(reference::A);
+ drawCircle(reference::B);
+ drawCircle(reference::C);
+ drawCircle(reference::Y);
+ drawCircle(reference::X);
+ drawCircle(reference::Z);
+ drawCircle(reference::U);
+ drawCircle(reference::K);
+ drawCircle(reference::M);
+ drawCircle(reference::G);
+ drawCircle(reference::H);
+ drawCircle(reference::J);
+ drawCircle(reference::I);
+ cr->restore();
+ }
+
+ // labels
+ {
+ auto const label = [&](LocalVector2D const &v, std::string l) {
+ moveTo(v);
+ cr->rel_move_to(0.005, -0.02);
+ cr->show_text(l);
+ };
+ auto font = Cairo::ToyFontFace::create(
+ "monospace", Cairo::FONT_SLANT_NORMAL, Cairo::FONT_WEIGHT_NORMAL);
+
+ cr->save();
+ cr->set_font_face(font);
+ cr->set_font_size(0.03);
+
+ label(A, "A");
+ label(B, "B");
+ label(C, "C");
+ label(D, "D");
+ label(X, "X");
+ label(Y, "Y");
+
+ cr->restore();
+ }
+
+ surface->write_to_png(filename);
+#endif
+}
+*/
+
+#include "trianglegeometry_tmpl.cc"
diff --git a/dune/tectonic/problem-data/grid/trianglegeometry.hh b/dune/tectonic/problem-data/grid/trianglegeometry.hh
new file mode 100644
index 00000000..92f178a6
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/trianglegeometry.hh
@@ -0,0 +1,88 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGEOMETRY_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGEOMETRY_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+// C
+// / |
+// / |
+// / |height
+// / |
+// A --- B
+// length
+
+template <class ctype = double>
+class TriangleGeometry {
+public:
+ typedef Dune::FieldVector<ctype, MY_DIM> GlobalCoords;
+ using WeakeningRegion = ConvexPolyhedron<GlobalCoords>;
+
+ static constexpr double lengthScale() {
+ return 1.0;
+ } // scaling factor
+
+private:
+ const ctype length_;
+ const ctype height_;
+#if MY_DIM == 3
+ const ctype depth_;
+#endif
+
+ // corners of triangle
+ const GlobalCoords A_;
+ const GlobalCoords B_;
+ const GlobalCoords C_;
+
+ // weakening regions
+ std::vector<WeakeningRegion> weakeningRegions_;
+
+public:
+#if MY_DIM == 3
+ TriangleGeometry(const GlobalCoords& origin,
+ const double length = 0.5, const double height = 0.5, const double depth = 0.12);
+
+ const auto& depth() const {
+ return depth_;
+ }
+#else
+ TriangleGeometry(const GlobalCoords& origin,
+ const double length = 0.5, const double height = 0.5);
+#endif
+
+ void addWeakeningRegion(const WeakeningRegion& weakeningRegion);
+ void addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1);
+
+ const auto& A() const {
+ return A_;
+ }
+
+ const auto& B() const {
+ return B_;
+ }
+
+ const auto& C() const {
+ return C_;
+ }
+
+ const auto& weakeningRegions() const {
+ return weakeningRegions_;
+ }
+
+ const auto& length() const {
+ return length_;
+ }
+
+ const auto& height() const {
+ return height_;
+ }
+
+ void write() const;
+
+ // void render() const;
+};
+#endif
diff --git a/dune/tectonic/problem-data/grid/trianglegeometry_tmpl.cc b/dune/tectonic/problem-data/grid/trianglegeometry_tmpl.cc
new file mode 100644
index 00000000..c67a450f
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/trianglegeometry_tmpl.cc
@@ -0,0 +1,8 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "trianglegeometry.hh"
+
+template class TriangleGeometry<typename Grid::ctype>;
diff --git a/dune/tectonic/problem-data/grid/trianglegrids.cc b/dune/tectonic/problem-data/grid/trianglegrids.cc
new file mode 100644
index 00000000..f45ada06
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/trianglegrids.cc
@@ -0,0 +1,204 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "trianglegrids.hh"
+#include "../midpoint.hh"
+#include "../../utils/diameter.hh"
+
+#include "simplexmanager.hh"
+
+// Fix: 3D case (still Elias code)
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<TriangleGeometry<typename Grid::ctype>>>& triangleGeometries) :
+ triangleGeometries_(triangleGeometries)
+{
+ size_t const gridCount = triangleGeometries.size();
+ grids.resize(gridCount);
+ gridFactories.resize(gridCount);
+
+ for (size_t idx=0; idx<grids.size(); idx++) {
+ const auto& triangleGeometry = *triangleGeometries[idx];
+
+ const auto& A = triangleGeometry.A();
+ const auto& B = triangleGeometry.B();
+ const auto& C = triangleGeometry.C();
+
+ const size_t vc = 3;
+
+#if MY_DIM == 3
+ Dune::FieldMatrix<double, 2 * vc, MY_DIM> vertices;
+#else
+ Dune::FieldMatrix<double, vc, MY_DIM> vertices;
+#endif
+ for (size_t i = 0; i < 2; ++i) {
+#if MY_DIM == 3
+ size_t numXYplanes = 2;
+#else
+ size_t numXYplanes = 1;
+#endif
+ size_t k = 0;
+ for (size_t j = 1; j <= numXYplanes; ++j) {
+ vertices[k++][i] = A[i];
+ vertices[k++][i] = B[i];
+ vertices[k++][i] = C[i];
+ assert(k == j * vc);
+ }
+ }
+
+#if MY_DIM == 3
+ for (size_t k = 0; k < vc; ++k) {
+ vertices[k][2] = 0;
+ vertices[k + vc][2] = triangleGeometry.depth();
+ }
+#endif
+
+ for (size_t i = 0; i < vertices.N(); ++i)
+ gridFactories[idx].insertVertex(vertices[i]);
+
+ Dune::GeometryType cell;
+#if MY_DIM == 3
+ cell = Dune::GeometryTypes::tetrahedron;
+#else
+ cell = Dune::GeometryTypes::triangle;
+#endif
+
+#if MY_DIM == 3
+ SimplexManager sm(vc);
+#else
+ SimplexManager sm;
+#endif
+ sm.addFromVerticesFFB(1, 2, 0);
+ auto const &simplices = sm.getSimplices();
+
+ // sanity-check choices of simplices
+ for (size_t i = 0; i < simplices.size(); ++i) {
+ Dune::FieldMatrix<double, MY_DIM, MY_DIM> check;
+ for (size_t j = 0; j < MY_DIM; ++j)
+ check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+ assert(check.determinant() > 0);
+ gridFactories[idx].insertElement(cell, simplices[i]);
+ }
+
+ grids[idx] = std::shared_ptr<Grid>(gridFactories[idx].createGrid());
+ grids[idx]->setRefinementType(Grid::RefinementType::COPY);
+ }
+}
+
+template <class Grid>
+std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() {
+ return grids;
+}
+
+template <class Grid>
+template <class GridView>
+MyFaces<GridView> GridsConstructor<Grid>::constructFaces(
+ GridView const &gridView, const TriangleGeometry<typename Grid::ctype>& triangleGeometry) {
+ return MyFaces<GridView>(gridView, triangleGeometry);
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+ Vector const &c) {
+ return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+ Vector const &x) {
+ auto const minmax0 = std::minmax(v1[0], v2[0]);
+ auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+ if (minmax0.first - 1e-14 * triangleGeometry_.lengthScale() > x[0] or
+ x[0] > minmax0.second + 1e-14 * triangleGeometry_.lengthScale())
+ return false;
+ if (minmax1.first - 1e-14 * triangleGeometry_.lengthScale() > x[1] or
+ x[1] > minmax1.second + 1e-14 * triangleGeometry_.lengthScale())
+ return false;
+
+ return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+ Vector const &x) {
+ return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+MyFaces<GridView>::MyFaces(GridView const &gridView, const TriangleGeometry<typename GridView::ctype>& triangleGeometry)
+ :
+ #if MY_DIM == 3
+ a(gridView),
+ b(gridView),
+ c(gridView),
+ top(gridView),
+ bottom(gridView),
+ #else
+ a(gridView),
+ b(gridView),
+ c(gridView),
+ #endif
+ triangleGeometry_(triangleGeometry)
+ {
+ a.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(triangleGeometry_.B(), triangleGeometry_.C(), in.geometry().center());
+ });
+
+ b.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(triangleGeometry_.A(), triangleGeometry_.C(), in.geometry().center());
+ });
+
+ c.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return xyBetween(triangleGeometry_.A(), triangleGeometry_.B(), in.geometry().center());
+ });
+
+ #if MY_DIM == 3
+ top.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return isClose(triangleGeometry_.depth(), in.geometry().center()[2]);
+ });
+
+ bottom.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return isClose(0.0, in.geometry().center()[2]);
+ });
+ #endif
+ }
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
+ return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
+}
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+ double smallestDiameter, double lengthScale) {
+ bool needRefine = true;
+ while (true) {
+ needRefine = false;
+ for (auto &&e : elements(grid.leafGridView())) {
+ auto const geometry = e.geometry();
+
+ auto const weakeningRegionDistance =
+ distance(weakPatch, geometry, 1e-6 * lengthScale);
+ auto const admissibleDiameter =
+ computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter, lengthScale);
+
+ if (diameter(geometry) <= admissibleDiameter)
+ continue;
+
+ needRefine = true;
+ grid.mark(1, e);
+ }
+ if (!needRefine)
+ break;
+
+ grid.preAdapt();
+ grid.adapt();
+ grid.postAdapt();
+ }
+}
+
+#include "trianglegrids_tmpl.cc"
diff --git a/dune/tectonic/problem-data/grid/trianglegrids.hh b/dune/tectonic/problem-data/grid/trianglegrids.hh
new file mode 100644
index 00000000..032cfa7d
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/trianglegrids.hh
@@ -0,0 +1,77 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGRIDS_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGRIDS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include "trianglegeometry.hh"
+
+// C
+// / |
+// / |
+// b / | a
+// / |
+// A --- B
+// c
+
+template <class GridView> struct MyFaces {
+ BoundaryPatch<GridView> a;
+ BoundaryPatch<GridView> b;
+ BoundaryPatch<GridView> c;
+
+#if MY_DIM == 3
+ BoundaryPatch<GridView> top;
+ BoundaryPatch<GridView> bottom;
+#endif
+
+ MyFaces(GridView const &gridView, const TriangleGeometry<typename GridView::ctype>& triangleGeometry);
+
+private:
+ const TriangleGeometry<typename GridView::ctype>& triangleGeometry_;
+
+ bool isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14 * triangleGeometry_.lengthScale();
+ }
+
+ bool isClose2(double a, double b) {
+ return std::abs(a - b) <
+ 1e-14 * triangleGeometry_.lengthScale() * triangleGeometry_.lengthScale();
+ }
+
+ template <class Vector>
+ bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
+
+ template <class Vector>
+ bool xyCollinear(Vector const &a, Vector const &b, Vector const &c);
+
+ template <class Vector>
+ bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
+};
+
+template <class Grid> class GridsConstructor {
+public:
+ GridsConstructor(const std::vector<std::shared_ptr<TriangleGeometry<typename Grid::ctype>>>& triangleGeometries);
+
+ std::vector<std::shared_ptr<Grid>>& getGrids();
+
+ template <class GridView>
+ MyFaces<GridView> constructFaces(const GridView& gridView, const TriangleGeometry<typename Grid::ctype>& triangleGeometry);
+
+private:
+ const std::vector<std::shared_ptr<TriangleGeometry<typename Grid::ctype>>>& triangleGeometries_;
+ std::vector<Dune::GridFactory<Grid>> gridFactories;
+ std::vector<std::shared_ptr<Grid>> grids;
+};
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale);
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+ double smallestDiameter, double lengthScale);
+
+#endif
+
+
diff --git a/dune/tectonic/problem-data/grid/trianglegrids_tmpl.cc b/dune/tectonic/problem-data/grid/trianglegrids_tmpl.cc
new file mode 100644
index 00000000..d76f933c
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/trianglegrids_tmpl.cc
@@ -0,0 +1,24 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+#include "trianglegeometry.hh"
+
+template class GridsConstructor<Grid>;
+
+template struct MyFaces<DeformedGrid::LeafGridView>;
+template struct MyFaces<DeformedGrid::LevelGridView>;
+
+template MyFaces<DeformedGrid::LeafGridView> GridsConstructor<Grid>::constructFaces(
+ const DeformedGrid::LeafGridView&,
+ const TriangleGeometry<typename Grid::ctype>&);
+
+template MyFaces<DeformedGrid::LevelGridView> GridsConstructor<Grid>::constructFaces(
+ const DeformedGrid::LevelGridView&,
+ const TriangleGeometry<typename Grid::ctype>&);
+
+template void refine<Grid, LocalVector>(
+ Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+ double smallestDiameter, double lengthScale);
diff --git a/dune/tectonic/problem-data/segmented-function.hh b/dune/tectonic/problem-data/segmented-function.hh
index 866415af..a3bc7b7d 100644
--- a/dune/tectonic/problem-data/segmented-function.hh
+++ b/dune/tectonic/problem-data/segmented-function.hh
@@ -15,10 +15,10 @@ class SegmentedFunction
Dune::FieldVector<double, MY_DIM> const &y,
Dune::FieldVector<double, MY_DIM> const &z) const {
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 false; //TODO liesBelow(MyGeometry::K, MyGeometry::M, z);
- };
+ }
double const _v1;
double const _v2;
diff --git a/dune/tectonic/problem-data/strikeslipbody.hh b/dune/tectonic/problem-data/strikeslipbody.hh
new file mode 100644
index 00000000..2e63a6af
--- /dev/null
+++ b/dune/tectonic/problem-data/strikeslipbody.hh
@@ -0,0 +1,89 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYBODYDATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include "../data-structures/body/bodydata.hh"
+#include "gravity.hh"
+#include "grid/cuboidgeometry.hh"
+
+
+class SegmentedFunction
+ : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+ Dune::FieldVector<double, 1>> {
+private:
+ bool distFromDiagonal(const Dune::FieldVector<double, MY_DIM>& z) const {
+ return std::abs(z[0] - z[1])/std::sqrt(2);
+ }
+
+ bool insideViscousRegion(Dune::FieldVector<double, MY_DIM> const &z) const {
+ return distFromDiagonal(z) > distFromDiag_;
+ }
+
+ const double distFromDiag_;
+
+ const double v1_;
+ const double v2_;
+
+public:
+ SegmentedFunction(double v1, double v2, double distFromDiag = 0.1) :
+ distFromDiag_(distFromDiag) , v1_(v1), v2_(v2) {}
+
+ void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+ Dune::FieldVector<double, 1> &y) const {
+ y = insideViscousRegion(x) ? v2_ : v1_;
+ }
+};
+
+
+template <int dimension> class MyBodyData : public BodyData<dimension> {
+ using typename BodyData<dimension>::ScalarFunction;
+ using typename BodyData<dimension>::VectorField;
+
+public:
+ MyBodyData(Dune::ParameterTree const &parset, const double gravity, const Dune::FieldVector<double, dimension>& zenith)
+ : poissonRatio_(parset.get<double>("poissonRatio")),
+ youngModulus_(3.0 * parset.get<double>("bulkModulus") *
+ (1.0 - 2.0 * poissonRatio_)),
+ zenith_(zenith),
+ shearViscosityField_(
+ parset.get<double>("elastic.shearViscosity"),
+ parset.get<double>("viscoelastic.shearViscosity"),
+ parset.get<double>("elastic.distFromDiag")),
+ bulkViscosityField_(
+ parset.get<double>("elastic.bulkViscosity"),
+ parset.get<double>("viscoelastic.bulkViscosity"),
+ parset.get<double>("elastic.distFromDiag")),
+ densityField_(parset.get<double>("elastic.density"),
+ parset.get<double>("viscoelastic.density"),
+ parset.get<double>("elastic.distFromDiag")),
+ gravityField_(densityField_, zenith_,
+ gravity) {}
+
+ double getPoissonRatio() const override { return poissonRatio_; }
+ double getYoungModulus() const override { return youngModulus_; }
+ ScalarFunction const &getShearViscosityField() const override {
+ return shearViscosityField_;
+ }
+ ScalarFunction const &getBulkViscosityField() const override {
+ return bulkViscosityField_;
+ }
+ ScalarFunction const &getDensityField() const override {
+ return densityField_;
+ }
+ VectorField const &getGravityField() const override { return gravityField_; }
+
+private:
+ double const poissonRatio_;
+ double const youngModulus_;
+ Dune::FieldVector<double, dimension> zenith_;
+
+ SegmentedFunction const shearViscosityField_;
+ SegmentedFunction const bulkViscosityField_;
+ SegmentedFunction const densityField_;
+ Gravity<dimension> const gravityField_;
+};
+#endif
diff --git a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
index 303e1ea0..c077470e 100644
--- a/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
+++ b/dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc
@@ -121,7 +121,7 @@ template <class field_type, class GridType0, class GridType1>
void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagrangeMatrix()
{
// clear matrix
- const auto nNonmortarVertices = nonmortarBoundary_.numVertices() - crosspoints_.size();
+ const auto nNonmortarVertices = nonmortarBoundary_.numVertices(); // - crosspoints_.size();
nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nNonmortarVertices);
nonmortarLagrangeMatrix_ = 0;
@@ -184,7 +184,7 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
hasObstacle_.resize(vertices.size());
int idx = 0;
for (size_t i=0; i<globalToLocal_.size(); i++) {
- bool nonmortarVertex = vertices[i][0] and crosspoints_.count(i)==0;
+ bool nonmortarVertex = vertices[i][0]; // and crosspoints_.count(i)==0;
globalToLocal_[i] = nonmortarVertex ? idx++ : -1;
hasObstacle_[i][0] = nonmortarVertex ? true : false;
}
@@ -193,7 +193,7 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
assembleNonmortarLagrangeMatrix();
// The weak obstacle vector
- const auto nNonmortarVertices = nonmortarBoundary_.numVertices()-crosspoints_.size();
+ const auto nNonmortarVertices = nonmortarBoundary_.numVertices();// - crosspoints_.size();
weakObstacle_.resize(nNonmortarVertices);
weakObstacle_ = 0;
@@ -254,7 +254,7 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
std::vector<Dune::FieldVector<ctype,dim> > avNormals;
avNormals = nonmortarBoundary_.getNormals();
- // loop over all intersections and compute the matrix entries
+ // loop over all intersections and compute the matrix entriescrosspoints_
for (const auto& rIs : intersections(*glue_)) {
const auto& inside = rIs.inside();
@@ -301,13 +301,15 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
int globalIdx = indexSet0.subIndex(inside,faceIdxi,dim);
- if (crosspoints_.count(globalIdx)>0) {
+ /*if (crosspoints_.count(globalIdx)>0) {
isCrosspointFace = true;
- } else {
+ } else {*/
nonmortarFaceNodes.push_back(faceIdxi);
- }
+ //}
}
+
+
for (const auto& quadPt : quadRule) {
// compute integration element of overlap
@@ -328,7 +330,7 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
//nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
- if (isCrosspointFace) {
+ /*if (isCrosspointFace) {
constantDualFE.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);
// loop over all Lagrange multiplier shape functions
@@ -355,7 +357,7 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
}
}
- } else {
+ } else {*/
dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
// loop over all Lagrange multiplier shape functions
@@ -382,7 +384,7 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
}
}
- }
+ //}
}
}
diff --git a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
index 0c0d7ca9..75ede91c 100644
--- a/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
+++ b/dune/tectonic/spatial-solving/contact/nbodyassembler.cc
@@ -15,7 +15,7 @@
template <class GridType, class VectorType>
void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
- /* std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
+ /*std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
for (size_t i=0; i<nGrids(); i++) {
bodywiseBoundaries[i] = *dirichletVertices_[i];
}
@@ -57,9 +57,9 @@ void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
}
- for (size_t i=0; i<crosspoints_.size(); i++) {
+ /*for (size_t i=0; i<crosspoints_.size(); i++) {
print(crosspoints_[i], "crosspoints " + std::to_string(i));
- }
+ }*/
}
template <class GridType, class VectorType>
@@ -94,18 +94,18 @@ void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
contactCoupling_[i]->reducePatches();
}
- setCrosspoints();
+ //setCrosspoints();
- for (int i=0; i<nGrids(); i++) {
+ /*for (int i=0; i<nGrids(); i++) {
print(crosspoints_[i], "crosspoints " + std::to_string(i));
- }
+ }*/
for (size_t i=0; i<contactCoupling_.size(); i++) {
contactCoupling_[i]->setCrosspoints(crosspoints_[coupling_[i].gridIdx_[0]]);
contactCoupling_[i]->setup();
- print(contactCoupling_[i]->nonmortarLagrangeMatrix() , "nonmortarLagrangeMatrix " + std::to_string(i) + ": ");
- print(contactCoupling_[i]->mortarLagrangeMatrix() , "mortarLagrangeMatrix " + std::to_string(i) + ": ");
+ //print(contactCoupling_[i]->nonmortarLagrangeMatrix() , "nonmortarLagrangeMatrix " + std::to_string(i) + ": ");
+ //print(contactCoupling_[i]->mortarLagrangeMatrix() , "mortarLagrangeMatrix " + std::to_string(i) + ": ");
}
// setup Householder reflections
@@ -122,17 +122,35 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
// Compute local coordinate systems for the dofs with obstacles
// //////////////////////////////////////////////////////////////////
- std::vector<std::vector<MatrixBlock> > coordinateSystems(nCouplings());
+ // first canonical basis vector
+ using CoordinateType = Dune::FieldVector<field_type,dim>;
+ CoordinateType e0(0);
+ e0[0] = 1;
+
+ // local identity
+ Dune::ScaledIdentityMatrix<field_type,dim> id(1);
+
+ // bodywise local coordinate systems: initialise with identity
+ std::vector<std::vector<MatrixBlock> > coordinateSystems(nGrids());
+ for (size_t i=0; i<coordinateSystems.size(); i++) {
+ coordinateSystems[i].resize(grids_[i]->size(dim));
+ for (size_t j=0; j<coordinateSystems[i].size(); j++) {
+ coordinateSystems[i][j] = id;
+ }
+ }
+
+ std::vector<std::map<size_t, CoordinateType>> crosspointDirections(nGrids());
for (int i=0; i<nCouplings(); i++) {
+ auto& coupling = coupling_[i];
- double dist = coupling_[i].obsDistance_;
- auto projection = coupling_[i].projection();
+ double dist = coupling.obsDistance_;
+ auto projection = coupling.projection();
if (!projection)
DUNE_THROW(Dune::Exception, "You have to supply a contact projection for the " << i << "-th coupling!");
- std::vector<Dune::FieldVector<field_type,dim> > directions;
+ std::vector<CoordinateType> directions;
const auto& nonmortarBoundary = contactCoupling_[i]->nonmortarBoundary();
const auto& mortarBoundary = contactCoupling_[i]->mortarBoundary();
@@ -140,7 +158,31 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
projection->project(nonmortarBoundary, mortarBoundary,dist);
projection->getObstacleDirections(directions);
- this->computeLocalCoordinateSystems(nonmortarBoundary,coordinateSystems[i], directions);
+ const auto nonmortarIdx = coupling_[i].gridIdx_[0];
+ auto globalToLocal = nonmortarBoundary.makeGlobalToLocal();
+
+ for (size_t j=0; j<globalToLocal.size(); j++) {
+ if (globalToLocal[j] > -1) {
+ // There is an obstacle at this vertex --> the coordinate system
+ // will consist of the surface normal and tangential vectors
+
+ if (crosspoints_[nonmortarIdx].count(j)>0) {
+ crosspointDirections[nonmortarIdx][j] += directions[globalToLocal[j]];
+ } else {
+ this->householderReflection(e0, directions[globalToLocal[j]], coordinateSystems[nonmortarIdx][j]);
+ }
+ }
+ }
+ }
+
+ // local coordinate systems at crosspoints
+ for (size_t i=0; i<crosspoints_.size(); i++) {
+ for (auto c : crosspoints_[i]) {
+ auto& direction = crosspointDirections[i][c];
+ direction /= direction.two_norm();
+
+ this->householderReflection(e0, direction, coordinateSystems[i][c]);
+ }
}
// ///////////////////////////////////////////////////////////////
@@ -148,28 +190,17 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
// ///////////////////////////////////////////////////////////////
this->localCoordSystems_.resize(numDofs());
- // initialise with identity
- Dune::ScaledIdentityMatrix<field_type,dim> id(1);
- for (size_t i=0; i<this->localCoordSystems_.size(); i++)
- this->localCoordSystems_[i] = id;
-
- for (int j=0; j<nCouplings(); j++) {
- int grid0Idx = coupling_[j].gridIdx_[0];
-
- // compute offset
- int idx = 0;
- for (int k=0;k<grid0Idx;k++)
- idx += grids_[k]->size(dim);
-
- const auto& globalToLocal = contactCoupling_[j]->globalToLocal();
+ size_t offSet = 0;
+ for (size_t i=0; i<nGrids(); i++) {
+ for (std::size_t k=0; k<coordinateSystems[i].size(); k++) {
+ this->localCoordSystems_[offSet+k] = coordinateSystems[i][k];
+ }
- // if a body is non-mortar more than once, one has to be careful with not overwriting entries
- for (std::size_t k=0; k<coordinateSystems[j].size(); k++)
- if(globalToLocal[k] > -1 ) //or crosspoints.count(k) > 0)
- this->localCoordSystems_[idx+k] = coordinateSystems[j][k];
+ offSet += grids_[i]->size(dim);
}
//print(this->localCoordSystems_, "localCoordSystems:");
+ //std::cout << "done" << std::endl;
}
template <class GridType, class VectorType>
@@ -239,7 +270,7 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
for (const auto& v : vertices(leafView)) {
int vIdx = indexSet.index(v);
- if (globalToLocal[vIdx]<0)
+ if (globalToLocal[vIdx]<0) // and crosspoints_[grid0Idx].count(vIdx)<1)
continue;
// Set the obstacle
@@ -252,6 +283,12 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
case Dune::Contact::CouplingPairBase::STICK_SLIP:
totalObstacles_[idx+vIdx].lower(0) = 0;
totalObstacles_[idx+vIdx].upper(0) = 0;
+
+ /*if (crosspoints_[grid0Idx].count(vIdx) > 0) {
+ totalObstacles_[idx+vIdx].lower(1) = 0;
+ totalObstacles_[idx+vIdx].upper(1) = 0;
+ }*/
+
break;
case Dune::Contact::CouplingPairBase::GLUE:
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.cc b/dune/tectonic/spatial-solving/fixedpointiterator.cc
index 39610f48..0a9ad5c1 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.cc
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.cc
@@ -39,6 +39,8 @@
#include "tnnmg/functional.hh"
#include "tnnmg/zerononlinearity.hh"
+#include <dune/tectonic/utils/reductionfactors.hh>
+
#include "fixedpointiterator.hh"
void FixedPointIterationCounter::operator+=(
@@ -225,7 +227,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
- relativeVelocities(v_m, v_rel);
+ split(v_m, v_rel);
//print(v_m, "v_m: ");
@@ -292,25 +294,21 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
}*/
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::relativeVelocities(
+void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::split(
const Vector& v,
- std::vector<Vector>& v_rel) const {
+ std::vector<Vector>& splitV) const {
const size_t nBodies = nBodyAssembler_.nGrids();
- // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
size_t globalIdx = 0;
- v_rel.resize(nBodies);
+ splitV.resize(nBodies);
for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
- const auto& nonmortarBoundary = *bodywiseNonmortarBoundaries_[bodyIdx];
- auto& v_rel_ = v_rel[bodyIdx];
+ auto& splitV_ = splitV[bodyIdx];
- v_rel_.resize(nonmortarBoundary.size());
+ splitV_.resize(nBodyAssembler_.grids_[bodyIdx]->size(dims));
- for (size_t i=0; i<v_rel_.size(); i++) {
- if (toBool(nonmortarBoundary[i])) {
- v_rel_[i] = v[globalIdx];
- }
+ for (size_t i=0; i<splitV_.size(); i++) {
+ splitV_[i] = v[globalIdx];
globalIdx++;
}
}
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.hh b/dune/tectonic/spatial-solving/fixedpointiterator.hh
index ae0912ad..183f2f8d 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.hh
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.hh
@@ -47,7 +47,7 @@ class FixedPointIterator {
using BitVector = Dune::BitSetVector<1>;
private:
- void relativeVelocities(const Vector& v, std::vector<Vector>& v_rel) const;
+ void split(const Vector& v, std::vector<Vector>& v_rel) const;
public:
FixedPointIterator(const Dune::ParameterTree& parset,
diff --git a/dune/tectonic/spatial-solving/tnnmg/linearization.hh b/dune/tectonic/spatial-solving/tnnmg/linearization.hh
index 556987c1..8abd074b 100644
--- a/dune/tectonic/spatial-solving/tnnmg/linearization.hh
+++ b/dune/tectonic/spatial-solving/tnnmg/linearization.hh
@@ -160,7 +160,7 @@ class Linearization {
//print(hessian_, "Hessian: ");
//std::cout << "--------" << (double) hessian_[21][21] << "------------" << std::endl;
- auto B = hessian_;
+ //auto B = hessian_;
//print(x, "x: ");
@@ -168,7 +168,7 @@ class Linearization {
phi.addHessian(x, hessian_);
//truncateMatrix(hessian_, regularityTruncationFlags, regularityTruncationFlags);
- B -= hessian_;
+ //B -= hessian_;
//print(B, "phi hessian: ");
// compute gradient
diff --git a/dune/tectonic/tests/CMakeLists.txt b/dune/tectonic/tests/CMakeLists.txt
index ef87c47f..69e6b8f1 100644
--- a/dune/tectonic/tests/CMakeLists.txt
+++ b/dune/tectonic/tests/CMakeLists.txt
@@ -1,3 +1,6 @@
+add_subdirectory("contacttest")
+add_subdirectory("hdf5test")
+
dune_add_test(SOURCES globalfrictioncontainertest.cc)
dune_add_test(SOURCES gridgluefrictiontest.cc)
dune_add_test(SOURCES nodalweightstest.cc)
diff --git a/dune/tectonic/tests/contacttest/CMakeLists.txt b/dune/tectonic/tests/contacttest/CMakeLists.txt
new file mode 100644
index 00000000..fc6644b8
--- /dev/null
+++ b/dune/tectonic/tests/contacttest/CMakeLists.txt
@@ -0,0 +1,33 @@
+add_custom_target(tectonic_tests_contacttest SOURCES
+ staticcontacttest.cfg
+ staticcontacttest-2D.cfg
+ staticcontacttest-3D.cfg
+)
+
+set(CONTACTTEST_SOURCE_FILES
+ ../../assemblers.cc
+ ../../data-structures/body/body.cc
+ ../../data-structures/network/levelcontactnetwork.cc
+ ../../data-structures/network/contactnetwork.cc
+ ../../data-structures/enumparser.cc
+ ../../factories/twoblocksfactory.cc
+ ../../factories/threeblocksfactory.cc
+ ../../factories/stackedblocksfactory.cc
+ ../../io/vtk.cc
+ ../../problem-data/grid/cuboidgeometry.cc
+ ../../problem-data/grid/mygrids.cc
+ ../../problem-data/grid/simplexmanager.cc
+ ../../spatial-solving/solverfactory.cc
+ staticcontacttest.cc
+)
+
+foreach(_dim 2 3)
+ set(_contacttest_target staticcontacttest-${_dim}D)
+
+ dune_add_test(NAME ${_contacttest_target} SOURCES ${CONTACTTEST_SOURCE_FILES})
+
+ add_dune_ug_flags(${_contacttest_target})
+ add_dune_hdf5_flags(${_contacttest_target})
+
+ set_property(TARGET ${_contacttest_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+endforeach()
diff --git a/dune/tectonic/tests/contacttest/staticcontacttest-2D.cfg b/dune/tectonic/tests/contacttest/staticcontacttest-2D.cfg
new file mode 100644
index 00000000..81b43367
--- /dev/null
+++ b/dune/tectonic/tests/contacttest/staticcontacttest-2D.cfg
@@ -0,0 +1,12 @@
+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter = 0.3 # 2e-3 [m]
+
+[u0.solver]
+tolerance = 1e-8
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance = 1e-10
diff --git a/dune/tectonic/tests/contacttest/staticcontacttest-3D.cfg b/dune/tectonic/tests/contacttest/staticcontacttest-3D.cfg
new file mode 100644
index 00000000..2e0af675
--- /dev/null
+++ b/dune/tectonic/tests/contacttest/staticcontacttest-3D.cfg
@@ -0,0 +1,15 @@
+# -*- mode:conf -*-
+[body0]
+smallestDiameter = 0.005 # 2e-3 [m]
+
+[body1]
+smallestDiameter = 0.005 # 2e-3 [m]
+
+[u0.solver]
+tolerance = 1e-8
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance = 1e-10
diff --git a/dune/tectonic/tests/contacttest/staticcontacttest.cc b/dune/tectonic/tests/contacttest/staticcontacttest.cc
new file mode 100644
index 00000000..c074fe1a
--- /dev/null
+++ b/dune/tectonic/tests/contacttest/staticcontacttest.cc
@@ -0,0 +1,313 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#define MY_DIM 2
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+
+#include <dune/fufem/formatstring.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/matrices.hh>
+
+#include <dune/tectonic/factories/twoblocksfactory.hh>
+#include <dune/tectonic/factories/threeblocksfactory.hh>
+#include <dune/tectonic/factories/stackedblocksfactory.hh>
+
+#include <dune/tectonic/io/vtk.hh>
+
+#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/zerononlinearity.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
+#include <dune/tectonic/spatial-solving/solverfactory.hh>
+#include <dune/tectonic/spatial-solving/contact/nbodycontacttransfer.hh>
+
+#include <dune/tectonic/utils/debugutils.hh>
+
+// for getcwd
+#include <unistd.h>
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/contacttest/staticcontacttest.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/contacttest/staticcontacttest-%dD.cfg", dims), parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+ return parset;
+}
+
+int main(int argc, char *argv[]) {
+ try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ char buffer[256];
+ char *val = getcwd(buffer, sizeof(buffer));
+ if (val) {
+ std::cout << buffer << std::endl;
+ std::cout << argv[0] << std::endl;
+ }
+
+ std::ofstream out("staticcontacttest.log");
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+ auto const parset = getParameters(argc, argv);
+
+ // ----------------------
+ // set up contact network
+ // ----------------------
+ using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
+ BlocksFactory blocksFactory(parset);
+
+ using ContactNetwork = typename BlocksFactory::ContactNetwork;
+ blocksFactory.build();
+
+ auto& contactNetwork = blocksFactory.contactNetwork();
+
+ const size_t bodyCount = contactNetwork.nBodies();
+
+ /* for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+ // printDofLocation(contactNetwork.body(i)->gridView());
+
+ //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+ //def = 1;
+ //deformedGridComplex.setDeformation(def, i);
+
+ const auto& level = *contactNetwork.level(i);
+
+ for (size_t j=0; j<level.nBodies(); j++) {
+ writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+ }
+ }*/
+
+ for (size_t i=0; i<bodyCount; i++) {
+ writeToVTK(contactNetwork.body(i)->gridView(), "", "initial_body_" + std::to_string(i) + "_leaf");
+ }
+
+ // ----------------------------
+ // assemble contactNetwork
+ // ----------------------------
+ contactNetwork.assemble();
+
+ //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+
+ // ----------------------------
+ // compute minimal stress solution
+ // ----------------------------
+ using Matrix = typename ContactNetwork::Matrix;
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+ // 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
+ BitVector dirichletNodes;
+ contactNetwork.totalNodes("dirichlet", dirichletNodes);
+
+ // minimal stress displacement
+ std::vector<Vector> u(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ u[i].resize(contactNetwork.body(i)->nVertices());
+ u[i] = 0.0;
+ }
+
+ // set rhs
+ std::vector<Vector> ell0(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ ell0[i].resize(u[i].size());
+ ell0[i] = 0.0;
+
+ contactNetwork.body(i)->externalForce()(0.0, ell0[i]);
+ }
+
+ // set elasticity operator
+ const auto& matrices = contactNetwork.matrices().elasticity;
+ std::vector<const Matrix*> matrices_ptr(matrices.size());
+ for (size_t i=0; i<matrices_ptr.size(); i++) {
+ matrices_ptr[i] = matrices[i].get();
+ }
+
+ // assemble full global contact problem
+ Matrix bilinearForm;
+
+ nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs;
+ nBodyAssembler.assembleRightHandSide(ell0, totalRhs);
+
+ Vector totalX;
+ nBodyAssembler.nodalToTransformed(u, totalX);
+
+ // get lower and upper obstacles
+ const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+ Vector lower(totalObstacles.size());
+ Vector upper(totalObstacles.size());
+
+ for (size_t j=0; j<totalObstacles.size(); ++j) {
+ const auto& totalObstaclesj = totalObstacles[j];
+ auto& lowerj = lower[j];
+ auto& upperj = upper[j];
+ for (size_t d=0; d<dims; ++d) {
+ lowerj[d] = totalObstaclesj[d][0];
+ upperj[d] = totalObstaclesj[d][1];
+ }
+ }
+
+ // print problem
+ print(bilinearForm, "bilinearForm");
+ print(totalRhs, "totalRhs");
+ print(dirichletNodes, "ignore");
+
+ std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+ contactNetwork.frictionNodes(frictionNodes);
+
+ for (size_t i=0; i<frictionNodes.size(); i++) {
+ print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+ }
+
+ print(totalObstacles, "totalObstacles");
+ print(lower, "lower");
+ print(upper, "upper");
+
+ // set up functional
+ using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+ Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
+
+ /* std::vector<BitVector> bodyDirichletNodes;
+ nBodyAssembler.postprocess(dirichletNodes, bodyDirichletNodes);
+ for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+ print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+ }*/
+
+ /* print(bilinearForm, "matrix: ");
+ print(totalX, "totalX: ");
+ print(totalRhs, "totalRhs: ");*/
+
+ // make linear solver for linear correction in TNNMGStep
+ const auto& solverParset = parset.sub("u0.solver");
+
+ using Norm = EnergyNorm<Matrix, Vector>;
+ using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+ // set multigrid solver
+ auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+ using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+ using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+ TransferOperators transfer(contactNetwork.nLevels()-1);
+ for (size_t i=0; i<transfer.size(); i++) {
+ transfer[i] = std::make_shared<TransferOperator>();
+ transfer[i]->setup(contactNetwork, i, i+1);
+ }
+
+ // Remove any recompute filed so that initially the full transferoperator is assembled
+ for (size_t i=0; i<transfer.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+ auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+ linearMultigridStep->setMGType(1, 3, 3);
+ linearMultigridStep->setSmoother(smoother);
+ linearMultigridStep->setTransferOperators(transfer);
+
+ Norm norm(*linearMultigridStep);
+
+ auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, parset.get<int>("solver.tnnmg.main.multi"), parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+
+ // set up TNMMG solver
+ using Factory = SolverFactory<Functional, BitVector>;
+ Factory factory(parset.sub("solver.tnnmg"), J, dirichletNodes);
+
+ factory.build(linearSolver);
+ auto tnnmgStep = factory.step();
+ factory.setProblem(totalX);
+
+ //const EnergyNorm<Matrix, Vector> norm(bilinearForm);
+
+ std::cout << "solving linear problem for u..." << std::endl;
+ LoopSolver<Vector> solver(
+ *tnnmgStep.get(), solverParset.get<size_t>("maximumIterations"),
+ solverParset.get<double>("tolerance"), norm,
+ solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+ solver.preprocess();
+ solver.solve();
+
+ nBodyAssembler.postprocess(tnnmgStep->getSol(), u);
+
+ print(tnnmgStep->getSol(), "totalX:");
+ print(u, "u:");
+
+ for (size_t i=0; i<bodyCount; i++) {
+ auto& body = contactNetwork.body(i);
+ body->setDeformation(u[i]);
+ writeToVTK(body->gridView(), "", "body_" + std::to_string(i) + "_leaf");
+ }
+
+ /*
+ using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+ using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+ using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+ StateUpdater<ScalarVector, Vector>>;
+
+ BoundaryFunctions velocityDirichletFunctions;
+ contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+ BoundaryNodes dirichletNodes;
+ contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+ for (size_t i=0; i<dirichletNodes.size(); i++) {
+ for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+ print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+ }
+ }*/
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+
+ } catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+ } catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+ }
+}
diff --git a/dune/tectonic/tests/contacttest/staticcontacttest.cfg b/dune/tectonic/tests/contacttest/staticcontacttest.cfg
new file mode 100644
index 00000000..f16fa5be
--- /dev/null
+++ b/dune/tectonic/tests/contacttest/staticcontacttest.cfg
@@ -0,0 +1,65 @@
+# -*- mode:conf -*-
+gravity = 9.81 # [m/s^2]
+
+[body]
+length = 0.4 # [m]
+height = 0.04 # [m]
+depth = 0.04 # [m]
+bulkModulus = 2190 # [Pa] #2190
+poissonRatio = 0.11 # [1] #0.11
+[body.elastic]
+density = 750 # [kg/m^3] #750
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+[body.viscoelastic]
+density = 1000 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+
+[boundary.friction]
+C = 6 # [Pa]
+mu0 = 0.48 # [ ]
+V0 = 1e-3 # [m/s]
+L = 1e-6 # [m]
+initialAlpha = 0 # [ ]
+stateModel = AgeingLaw
+frictionModel = Truncated #Regularised
+[boundary.friction.weakening]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
+[boundary.friction.strengthening]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
+
+[boundary.neumann]
+sigmaN = 0.0 # [Pa]
+
+[boundary.dirichlet]
+finalVelocity = 5e-5 # [m/s]
+
+[problem]
+bodyCount = 2
+
+[u0.solver]
+maximumIterations = 100
+verbosity = full
+
+[solver.tnnmg.preconditioner]
+mode = additive
+patchDepth = 1
+maximumIterations = 2
+verbosity = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity = quiet
+
+[solver.tnnmg.main]
+pre = 1
+multi = 5 # number of multigrid steps
+post = 0
+
+
+
diff --git a/dune/tectonic/tests/hdf5test/CMakeLists.txt b/dune/tectonic/tests/hdf5test/CMakeLists.txt
new file mode 100644
index 00000000..27378cc4
--- /dev/null
+++ b/dune/tectonic/tests/hdf5test/CMakeLists.txt
@@ -0,0 +1,33 @@
+add_custom_target(tectonic_tests_hdf5test SOURCES
+ hdf5test.cfg
+ hdf5test-2D.cfg
+)
+
+set(HDF5TEST_SOURCE_FILES
+ ../../assemblers.cc
+ ../../data-structures/body/body.cc
+ ../../data-structures/network/levelcontactnetwork.cc
+ ../../data-structures/network/contactnetwork.cc
+ ../../data-structures/enumparser.cc
+ ../../factories/strikeslipfactory.cc
+ ../../io/vtk.cc
+ ../../io/hdf5/frictionalboundary-writer.cc
+ ../../io/hdf5/iteration-writer.cc
+ ../../io/hdf5/time-writer.cc
+ ../../problem-data/grid/simplexmanager.cc
+ ../../problem-data/grid/trianglegeometry.cc
+ ../../problem-data/grid/trianglegrids.cc
+ ../../spatial-solving/solverfactory.cc
+ hdf5test.cc
+)
+
+foreach(_dim 2)
+ set(_hdf5test_target hdf5test-${_dim}D)
+
+ dune_add_test(NAME ${_hdf5test_target} SOURCES ${HDF5TEST_SOURCE_FILES})
+
+ add_dune_ug_flags(${_hdf5test_target})
+ add_dune_hdf5_flags(${_hdf5test_target})
+
+ set_property(TARGET ${_hdf5test_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+endforeach()
diff --git a/dune/tectonic/tests/hdf5test/hdf5test-2D.cfg b/dune/tectonic/tests/hdf5test/hdf5test-2D.cfg
new file mode 100644
index 00000000..f212cdea
--- /dev/null
+++ b/dune/tectonic/tests/hdf5test/hdf5test-2D.cfg
@@ -0,0 +1,27 @@
+# -*- mode:conf -*-
+[body0]
+smallestDiameter = 0.3 # 2e-3 [m]
+
+[body1]
+smallestDiameter = 0.3 # 2e-3 [m]
+
+[timeSteps]
+refinementTolerance = 1e-5 # 1e-5
+
+[u0.solver]
+tolerance = 1e-8
+
+[a0.solver]
+tolerance = 1e-8
+
+[v.solver]
+tolerance = 1e-8
+
+[v.fpi]
+tolerance = 1e-5
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance = 1e-10
diff --git a/dune/tectonic/tests/hdf5test/hdf5test.cc b/dune/tectonic/tests/hdf5test/hdf5test.cc
new file mode 100644
index 00000000..4953ae6e
--- /dev/null
+++ b/dune/tectonic/tests/hdf5test/hdf5test.cc
@@ -0,0 +1,202 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#define MY_DIM 2
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/program_state.hh>
+
+#include <dune/tectonic/factories/strikeslipfactory.hh>
+
+#include <dune/tectonic/io/vtk.hh>
+#include <dune/tectonic/io/hdf5-writer.hh>
+
+#include <dune/tectonic/utils/geocoordinate.hh>
+#include <dune/tectonic/utils/debugutils.hh>
+
+// for getcwd
+#include <unistd.h>
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/hdf5test/hdf5test.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/hdf5test/hdf5test-%dD.cfg", dims), parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+ return parset;
+}
+
+int main(int argc, char *argv[]) {
+ try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ char buffer[256];
+ char *val = getcwd(buffer, sizeof(buffer));
+ if (val) {
+ std::cout << buffer << std::endl;
+ std::cout << argv[0] << std::endl;
+ }
+
+ std::ofstream out("hdf5test.log");
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+ auto const parset = getParameters(argc, argv);
+
+ // ----------------------
+ // set up contact network
+ // ----------------------
+ using BlocksFactory = StrikeSlipFactory<Grid, Vector>;
+ BlocksFactory blocksFactory(parset);
+
+ using ContactNetwork = typename BlocksFactory::ContactNetwork;
+ blocksFactory.build();
+
+ auto& contactNetwork = blocksFactory.contactNetwork();
+
+ const size_t bodyCount = contactNetwork.nBodies();
+
+ /* for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+ // printDofLocation(contactNetwork.body(i)->gridView());
+
+ //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+ //def = 1;
+ //deformedGridComplex.setDeformation(def, i);
+
+ const auto& level = *contactNetwork.level(i);
+
+ for (size_t j=0; j<level.nBodies(); j++) {
+ writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+ }
+ }*/
+
+ for (size_t i=0; i<bodyCount; i++) {
+ writeToVTK(contactNetwork.body(i)->gridView(), "", "initial_body_" + std::to_string(i) + "_leaf");
+ }
+
+ // ----------------------------
+ // assemble contactNetwork
+ // ----------------------------
+ contactNetwork.assemble();
+
+ //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+ std::vector<size_t> nVertices(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ nVertices[i] = contactNetwork.body(i)->nVertices();
+ }
+
+ using MyProgramState = ProgramState<Vector, ScalarVector>;
+ MyProgramState programState(nVertices);
+ programState.setupInitialConditions(parset, contactNetwork);
+
+ auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ for (size_t i=0; i<bodyCount; i++) {
+ contactNetwork.body(i)->setDeformation(programState.u[i]);
+ }
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();
+
+ // ------------------------
+ // assemble global friction
+ // ------------------------
+ contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+
+ auto& globalFriction = contactNetwork.globalFriction();
+ globalFriction.updateAlpha(programState.alpha);
+
+ // -----------------
+ // init input/output
+ // -----------------
+ using Assembler = MyAssembler<DefLeafGridView, dims>;
+ using MyVertexBasis = typename Assembler::VertexBasis;
+ using MyCellBasis = typename Assembler::CellBasis;
+ std::vector<Vector> vertexCoordinates(bodyCount);
+ std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+ std::vector<const MyCellBasis* > cellBases(bodyCount);
+
+ auto& wPatches = blocksFactory.weakPatches();
+ std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ vertexBases[i] = &(body->assembler()->vertexBasis);
+ cellBases[i] = &(body->assembler()->cellBasis);
+
+ weakPatches[i].resize(1);
+ weakPatches[i][0] = wPatches[i].get();
+
+ auto& vertexCoords = vertexCoordinates[i];
+ vertexCoords.resize(nVertices[i]);
+
+ auto hostLeafView = body->grid()->hostGrid().leafGridView();
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(hostLeafView))
+ vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+ }
+
+ std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
+ contactNetwork.frictionPatches(frictionPatches);
+
+ auto const writeData = parset.get<bool>("io.data.write");
+ auto dataFile = writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
+ auto dataWriter =
+ writeData ? std::make_unique<
+ HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+ *dataFile, vertexCoordinates, vertexBases,
+ frictionPatches, weakPatches)
+ : nullptr;
+
+ IterationRegister iterationCount;
+
+ auto const report = [&](bool initial = false) {
+ if (writeData) {
+ dataWriter->reportSolution(programState, contactNetwork, globalFriction);
+ if (!initial)
+ dataWriter->reportIterations(programState, iterationCount);
+ dataFile->flush();
+ }
+ };
+ report(true);
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+
+ } catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+ } catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+ }
+}
diff --git a/dune/tectonic/tests/hdf5test/hdf5test.cfg b/dune/tectonic/tests/hdf5test/hdf5test.cfg
new file mode 100644
index 00000000..818e378d
--- /dev/null
+++ b/dune/tectonic/tests/hdf5test/hdf5test.cfg
@@ -0,0 +1,105 @@
+# -*- mode:conf -*-
+gravity = 9.81 # [m/s^2]
+
+[body0]
+length = 0.5 # [m]
+height = 0.5 # [m]
+depth = 0.12 # [m]
+bulkModulus = 1.5e5 # [Pa] #2190
+poissonRatio = 0.11 # [1] #0.11
+[body0.elastic]
+density = 1300 # [kg/m^3] #750
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+[body0.viscoelastic]
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+
+[body1]
+length = 0.5 # [m]
+height = 0.5 # [m]
+depth = 0.12 # [m]
+bulkModulus = 1.5e5 # [Pa]
+poissonRatio = 0.11 # [1]
+[body1.elastic]
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+[body1.viscoelastic]
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+
+[boundary.friction]
+C = 6 # [Pa]
+mu0 = 0.48 # [ ]
+V0 = 1e-3 # [m/s]
+L = 1e-6 # [m]
+initialAlpha = 0 # [ ]
+stateModel = AgeingLaw
+frictionModel = Truncated #Regularised
+[boundary.friction.weakening]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
+[boundary.friction.strengthening]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
+
+[boundary.neumann]
+sigmaN = 200.0 # [Pa]
+
+[boundary.dirichlet]
+finalVelocity = 1e-5 # [m/s]
+
+[io]
+data.write = true
+
+[problem]
+finalTime = 1000 # [s] #1000
+bodyCount = 2
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 2e-2 # 1e-6
+
+[timeSteps]
+scheme = newmark
+timeSteps = 1000
+
+[u0.solver]
+maximumIterations = 100
+verbosity = full
+
+[a0.solver]
+maximumIterations = 100
+verbosity = full
+
+[v.solver]
+maximumIterations = 100
+verbosity = quiet
+
+[v.fpi]
+maximumIterations = 10000
+lambda = 0.5
+
+[solver.tnnmg.preconditioner]
+mode = additive
+patchDepth = 1
+maximumIterations = 2
+verbosity = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity = quiet
+
+[solver.tnnmg.main]
+pre = 1
+multi = 5 # number of multigrid steps
+post = 0
+
+
+
diff --git a/dune/tectonic/time-stepping/adaptivetimestepper.cc b/dune/tectonic/time-stepping/adaptivetimestepper.cc
index 79bad14f..4c5b9185 100644
--- a/dune/tectonic/time-stepping/adaptivetimestepper.cc
+++ b/dune/tectonic/time-stepping/adaptivetimestepper.cc
@@ -8,8 +8,81 @@
#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+#include <dune/tectonic/utils/reductionfactors.hh>
+
#include "adaptivetimestepper.hh"
+template<class IterationStepType, class NormType, class ReductionFactorContainer>
+Dune::Solvers::Criterion reductionFactorCriterion(
+ IterationStepType& iterationStep,
+ const NormType& norm,
+ ReductionFactorContainer& reductionFactors)
+{
+ double normOfOldCorrection = 1;
+ auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+ return Dune::Solvers::Criterion(
+ [&, lastIterate, normOfOldCorrection] () mutable {
+ double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+ double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+
+ if (convRate>1.0)
+ std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+ normOfOldCorrection = normOfCorrection;
+ *lastIterate = *iterationStep.getIterate();
+
+ reductionFactors.push_back(convRate);
+ return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+ },
+ " reductionFactor");
+}
+
+
+template<class IterationStepType, class Functional, class ReductionFactorContainer>
+Dune::Solvers::Criterion energyCriterion(
+ const IterationStepType& iterationStep,
+ const Functional& f,
+ ReductionFactorContainer& reductionFactors)
+{
+ double normOfOldCorrection = 1;
+ auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+ return Dune::Solvers::Criterion(
+ [&, lastIterate, normOfOldCorrection] () mutable {
+ double normOfCorrection = std::abs(f(*lastIterate) - f(*iterationStep.getIterate())); //norm.diff(*lastIterate, *iterationStep.getIterate());
+
+ double convRate = (normOfOldCorrection != 0.0) ? 1.0 - (normOfCorrection / normOfOldCorrection) : 0.0;
+
+ if (convRate>1.0)
+ std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+ normOfOldCorrection = normOfCorrection;
+ *lastIterate = *iterationStep.getIterate();
+
+ reductionFactors.push_back(convRate);
+ return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+ },
+ " reductionFactor");
+}
+
+template <class ReductionFactorContainer>
+void updateReductionFactors(ReductionFactorContainer& reductionFactors) {
+ const size_t s = reductionFactors.size();
+
+ //print(reductionFactors, "reduction factors: ");
+
+ if (s>allReductionFactors.size()) {
+ allReductionFactors.resize(s);
+ }
+
+ for (size_t i=0; i<reductionFactors.size(); i++) {
+ allReductionFactors[i].push_back(reductionFactors[i]);
+ }
+
+ reductionFactors.clear();
+}
+
void IterationRegister::registerCount(FixedPointIterationCounter count) {
totalCount += count;
}
@@ -88,7 +161,7 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
of size tau/2 with one of size tau. The method makes multiple
coarsening/refining attempts, with coarsening coming first. */
- std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
+ //std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
// patch preconditioner only needs to be computed once per advance()
// make linear solver for linear correction in TNNMGStep
@@ -109,7 +182,8 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
Norm norm(*cgStep);
auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
-*/
+ */
+
// set multigrid solver
auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
@@ -135,13 +209,26 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+ Vector x;
+ x.resize(contactNetwork_.nBodyAssembler().totalHasObstacle_.size());
+ x = 0;
+ dynamic_cast<Dune::Solvers::IterationStep<Vector>*>(linearMultigridStep.get())->setProblem(x);
+ //dynamic_cast<Dune::Solvers::IterationStep<Vector>*>(cgStep.get())->setProblem(x);
+
const auto& currentNBodyAssembler = contactNetwork_.nBodyAssembler();
+ std::vector<double> reductionFactors;
+ linearSolver->addCriterion(reductionFactorCriterion(*linearMultigridStep, norm, reductionFactors));
+ //linearSolver->addCriterion(reductionFactorCriterion(*cgStep, norm, reductionFactors));
+
if (R1_.updaters == Updaters()) {
//setDeformation(current_);
R1_ = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, relativeTau_);
+ updateReductionFactors(reductionFactors);
}
+
+
//std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
size_t coarseningCount = 0;
@@ -157,6 +244,7 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
setDeformation(current_);
C = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, 2 * relativeTau_);
+ updateReductionFactors(reductionFactors);
std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
/*using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
@@ -167,6 +255,7 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
setDeformation(R1_.updaters);
auto&& nBodyAssembler = step(currentNBodyAssembler);
R2 = step(R1_.updaters, nBodyAssembler, linearSolver, relativeTime_ + relativeTau_, relativeTau_);
+ updateReductionFactors(reductionFactors);
std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
@@ -191,12 +280,14 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
while (true) {
setDeformation(current_);
F1 = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, relativeTau_ / 2.0);
+ updateReductionFactors(reductionFactors);
std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
setDeformation(F1.updaters);
auto&& nBodyAssembler = step(currentNBodyAssembler);
F2 = step(F1.updaters, nBodyAssembler, linearSolver, relativeTime_ + relativeTau_ / 2.0,
relativeTau_ / 2.0);
+ updateReductionFactors(reductionFactors);
std::cout << "AdaptiveTimeStepper F2 computed!" << std::endl << std::endl;
if (!mustRefine_(R1_.updaters, F2.updaters)) {
std::cout << "Sufficiently refined!" << std::endl;
diff --git a/dune/tectonic/time-stepping/state/ageinglawstateupdater.cc b/dune/tectonic/time-stepping/state/ageinglawstateupdater.cc
index 1a7dcea3..3d1feef2 100644
--- a/dune/tectonic/time-stepping/state/ageinglawstateupdater.cc
+++ b/dune/tectonic/time-stepping/state/ageinglawstateupdater.cc
@@ -54,10 +54,14 @@ auto liftSingularity(
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::solve(const Vector& velocity_field) {
- for (size_t i=0; i<alpha_.size(); i++) {
- double const V = velocity_field[i].two_norm();
- double const mtoL = -tau_ / L_;
- alpha_[i] = std::log(std::exp(alpha_o_[i] + V * mtoL) +
+
+ for (size_t i=0; i<alpha_.size(); ++i) {
+ auto tangentVelocity = velocity_field[localToGlobal_[i]];
+ tangentVelocity[0] = 0;
+
+ double const V = tangentVelocity.two_norm();
+ double const mtoL = -tau_ / L_;
+ alpha_[i] = std::log(std::exp(alpha_o_[i] + V * mtoL) +
V0_ * liftSingularity(mtoL, V));
}
}
diff --git a/dune/tectonic/time-stepping/state/sliplawstateupdater.cc b/dune/tectonic/time-stepping/state/sliplawstateupdater.cc
index 3b06f055..a024744c 100644
--- a/dune/tectonic/time-stepping/state/sliplawstateupdater.cc
+++ b/dune/tectonic/time-stepping/state/sliplawstateupdater.cc
@@ -41,8 +41,11 @@ template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::solve(
const Vector& velocity_field) {
- for (size_t i=0; i<alpha_.size(); i++) {
- double const V = velocity_field[i].two_norm();
+ for (size_t i=0; i<alpha_.size(); ++i) {
+ auto tangentVelocity = velocity_field[localToGlobal_[i]];
+ tangentVelocity[0] = 0;
+
+ double const V = tangentVelocity.two_norm();
double const mtVoL = -tau_ * V / L_;
alpha_[i] = (V <= 0) ? alpha_o_[i] : std::expm1(mtVoL) * std::log(V / V0_) +
alpha_o_[i] * std::exp(mtVoL);
diff --git a/dune/tectonic/utils/reductionfactors.hh b/dune/tectonic/utils/reductionfactors.hh
new file mode 100644
index 00000000..cc813928
--- /dev/null
+++ b/dune/tectonic/utils/reductionfactors.hh
@@ -0,0 +1 @@
+extern std::vector<std::vector<double>> allReductionFactors;
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e41c1f9a..080ed358 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,6 @@
add_subdirectory("foam")
add_subdirectory("multi-body-problem")
+add_subdirectory("strikeslip")
set(UGW_SOURCE_FILES
../dune/tectonic/assemblers.cc # FIXME
diff --git a/src/foam/foam-2D.cfg b/src/foam/foam-2D.cfg
index 1d8c5160..d911706c 100644
--- a/src/foam/foam-2D.cfg
+++ b/src/foam/foam-2D.cfg
@@ -1,9 +1,9 @@
# -*- mode:conf -*-
[body0]
-smallestDiameter = 0.02 # 2e-3 [m]
+smallestDiameter = 0.05 # 2e-3 [m]
[body1]
-smallestDiameter = 0.02 # 2e-3 [m]
+smallestDiameter = 0.01 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index 9e4c7e56..d78fc270 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -84,6 +84,9 @@
//#include <tbb/tbb.h> //TODO multi threading preconditioner?
//#include <pthread.h>
+#include <dune/tectonic/utils/reductionfactors.hh>
+std::vector<std::vector<double>> allReductionFactors;
+
size_t const dims = MY_DIM;
Dune::ParameterTree getParameters(int argc, char *argv[]) {
@@ -121,17 +124,14 @@ int main(int argc, char *argv[]) {
// ----------------------
// set up contact network
// ----------------------
- TwoBlocksFactory<Grid, Vector> twoBlocksFactory(parset);
- using ContactNetwork = typename TwoBlocksFactory<Grid, Vector>::ContactNetwork;
- twoBlocksFactory.build();
+ using BlocksFactory = TwoBlocksFactory<Grid, Vector>;
+ BlocksFactory blocksFactory(parset);
+ using ContactNetwork = typename BlocksFactory::ContactNetwork;
+ blocksFactory.build();
- ContactNetwork& contactNetwork = twoBlocksFactory.contactNetwork();
+ ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
- /*ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
- using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
- threeBlocksFactory.build();
- ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork(); */
const size_t bodyCount = contactNetwork.nBodies();
@@ -219,7 +219,7 @@ int main(int argc, char *argv[]) {
std::vector<const MyVertexBasis* > vertexBases(bodyCount);
std::vector<const MyCellBasis* > cellBases(bodyCount);
- auto& wPatches = twoBlocksFactory.weakPatches();
+ auto& wPatches = blocksFactory.weakPatches();
std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
@@ -234,20 +234,24 @@ int main(int argc, char *argv[]) {
auto& vertexCoords = vertexCoordinates[i];
vertexCoords.resize(nVertices[i]);
+ auto hostLeafView = body->grid()->hostGrid().leafGridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
- DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
- for (auto &&v : vertices(body->gridView()))
+ LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(hostLeafView))
vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
typename ContactNetwork::BoundaryPatches frictionBoundaries;
contactNetwork.boundaryPatches("friction", frictionBoundaries);
+ std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
+ contactNetwork.frictionPatches(frictionPatches);
+
auto dataWriter =
writeData ? std::make_unique<
HDF5LevelWriter<MyProgramState, MyVertexBasis, DefLeafGridView>>(
*dataFile, vertexCoordinates, vertexBases,
- frictionBoundaries, weakPatches)
+ frictionPatches, weakPatches)
: nullptr;
const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
@@ -256,7 +260,7 @@ int main(int argc, char *argv[]) {
auto const report = [&](bool initial = false) {
if (writeData) {
- dataWriter->reportSolution(programState, globalFriction);
+ dataWriter->reportSolution(programState, contactNetwork, globalFriction);
if (!initial)
dataWriter->reportIterations(programState, iterationCount);
dataFile->flush();
@@ -327,20 +331,20 @@ int main(int argc, char *argv[]) {
BoundaryNodes dirichletNodes;
contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- for (size_t i=0; i<dirichletNodes.size(); i++) {
+ /*for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
}
- }
+ }*/
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
contactNetwork.frictionNodes(frictionNodes);
- for (size_t i=0; i<frictionNodes.size(); i++) {
+ /*for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
- }
+ }*/
- DUNE_THROW(Dune::Exception, "Just need to stop here!");
+ //DUNE_THROW(Dune::Exception, "Just need to stop here!");
Updaters current(
initRateUpdater(
diff --git a/src/foam/foam.cfg b/src/foam/foam.cfg
index c8604d62..b5f4216c 100644
--- a/src/foam/foam.cfg
+++ b/src/foam/foam.cfg
@@ -2,52 +2,55 @@
gravity = 9.81 # [m/s^2]
[body0]
-length = 0.08 # [m]
+length = 0.4 # [m]
height = 0.04 # [m]
depth = 0.04 # [m]
-bulkModulus = 0.5e5 # [Pa] #2190
-poissonRatio = 0.30 # [1] #0.11
+bulkModulus = 1.5e5 # [Pa] #2190
+poissonRatio = 0.11 # [1] #0.11
[body0.elastic]
-density = 900 # [kg/m^3] #750
-shearViscosity = 1e20 # [Pas]
-bulkViscosity = 1e20 # [Pas]
+density = 1300 # [kg/m^3] #750
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
[body0.viscoelastic]
-density = 1000 # [kg/m^3]
-shearViscosity = 1e20 # [Pas]
-bulkViscosity = 1e20 # [Pas]
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
[body1]
length = 0.04 # [m]
height = 0.04 # [m]
depth = 0.04 # [m]
-bulkModulus = 0.5e5 # [Pa]
-poissonRatio = 0.30 # [1]
+bulkModulus = 1.5e5 # [Pa]
+poissonRatio = 0.11 # [1]
[body1.elastic]
-density = 900 # [kg/m^3]
-shearViscosity = 1e20 # [Pas]
-bulkViscosity = 1e20 # [Pas]
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
[body1.viscoelastic]
-density = 1000 # [kg/m^3]
-shearViscosity = 1e20 # [Pas]
-bulkViscosity = 1e20 # [Pas]
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
[boundary.friction]
-C = 10 # [Pa]
-mu0 = 0.7 # [ ]
-V0 = 5e-5 # [m/s]
-L = 2.25e-5 # [m]
+C = 6 # [Pa]
+mu0 = 0.48 # [ ]
+V0 = 1e-3 # [m/s]
+L = 1e-6 # [m]
initialAlpha = 0 # [ ]
stateModel = AgeingLaw
frictionModel = Truncated #Regularised
[boundary.friction.weakening]
-a = 0.002 # [ ]
-b = 0.017 # [ ]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
[boundary.friction.strengthening]
-a = 0.020 # [ ]
-b = 0.005 # [ ]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
[boundary.neumann]
-sigmaN = 0.0 # [Pa]
+sigmaN = 10000.0 # [Pa]
+
+[boundary.dirichlet]
+finalVelocity = 3e-3 # [m/s]
[io]
data.write = false
@@ -64,11 +67,11 @@ bodyCount = 2
[initialTime]
timeStep = 0
relativeTime = 0.0
-relativeTau = 5e-4 # 1e-6
+relativeTau = 5e-6 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 100000
+timeSteps = 1
[u0.solver]
maximumIterations = 100
diff --git a/src/multi-body-problem/multi-body-problem-2D.cfg b/src/multi-body-problem/multi-body-problem-2D.cfg
index cfee4140..641078a7 100644
--- a/src/multi-body-problem/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 0.1 # 2e-3 [m]
+smallestDiameter = 0.05 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem/multi-body-problem.cc b/src/multi-body-problem/multi-body-problem.cc
index eb2fd28f..62da56ba 100644
--- a/src/multi-body-problem/multi-body-problem.cc
+++ b/src/multi-body-problem/multi-body-problem.cc
@@ -87,6 +87,9 @@
size_t const dims = MY_DIM;
+#include <dune/tectonic/utils/reductionfactors.hh>
+std::vector<std::vector<double>> allReductionFactors;
+
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem/multi-body-problem.cfg", parset);
@@ -122,17 +125,15 @@ int main(int argc, char *argv[]) {
// ----------------------
// set up contact network
// ----------------------
- StackedBlocksFactory<Grid, Vector> stackedBlocksFactory(parset);
- using ContactNetwork = typename StackedBlocksFactory<Grid, Vector>::ContactNetwork;
- stackedBlocksFactory.build();
- ContactNetwork& contactNetwork = stackedBlocksFactory.contactNetwork();
+ //using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
+ using BlocksFactory = ThreeBlocksFactory<Grid, Vector>;
+ BlocksFactory blocksFactory(parset);
- /*ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
- using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
- threeBlocksFactory.build();
+ using ContactNetwork = typename BlocksFactory::ContactNetwork;
+ blocksFactory.build();
- ContactNetwork& contactNetwork = threeBlocksFactory.contactNetwork();*/
+ ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
const size_t bodyCount = contactNetwork.nBodies();
@@ -150,10 +151,10 @@ int main(int argc, char *argv[]) {
}*/
}
- for (size_t i=0; i<bodyCount; i++) {
+ /*for (size_t i=0; i<bodyCount; i++) {
printDofLocation(contactNetwork.body(i)->gridView());
writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
- }
+ } */
// ----------------------------
// assemble contactNetwork
@@ -223,8 +224,7 @@ int main(int argc, char *argv[]) {
std::vector<const MyVertexBasis* > vertexBases(bodyCount);
std::vector<const MyCellBasis* > cellBases(bodyCount);
- auto& wPatches = stackedBlocksFactory.weakPatches();
- //auto& wPatches = threeBlocksFactory.weakPatches();
+ auto& wPatches = blocksFactory.weakPatches();
std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
@@ -239,20 +239,24 @@ int main(int argc, char *argv[]) {
auto& vertexCoords = vertexCoordinates[i];
vertexCoords.resize(nVertices[i]);
+ auto hostLeafView = body->grid()->hostGrid().leafGridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
- DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
- for (auto &&v : vertices(body->gridView()))
+ LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(hostLeafView))
vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
typename ContactNetwork::BoundaryPatches frictionBoundaries;
contactNetwork.boundaryPatches("friction", frictionBoundaries);
+ std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
+ contactNetwork.frictionPatches(frictionPatches);
+
auto dataWriter =
writeData ? std::make_unique<
HDF5LevelWriter<MyProgramState, MyVertexBasis, DefLeafGridView>>(
*dataFile, vertexCoordinates, vertexBases,
- frictionBoundaries, weakPatches)
+ frictionPatches, weakPatches)
: nullptr;
const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
@@ -261,7 +265,7 @@ int main(int argc, char *argv[]) {
auto const report = [&](bool initial = false) {
if (writeData) {
- dataWriter->reportSolution(programState, globalFriction);
+ dataWriter->reportSolution(programState, contactNetwork, globalFriction);
if (!initial)
dataWriter->reportIterations(programState, iterationCount);
dataFile->flush();
@@ -317,7 +321,7 @@ int main(int argc, char *argv[]) {
}
}
- print(totalDirichletNodes, "totalDirichletNodes:");
+ //print(totalDirichletNodes, "totalDirichletNodes:");
//using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
@@ -334,18 +338,18 @@ int main(int argc, char *argv[]) {
BoundaryNodes dirichletNodes;
contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
- for (size_t i=0; i<dirichletNodes.size(); i++) {
+ /*for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
}
- }
+ }*/
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
contactNetwork.frictionNodes(frictionNodes);
- for (size_t i=0; i<frictionNodes.size(); i++) {
+ /*for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
- }
+ }*/
Updaters current(
initRateUpdater(
@@ -476,6 +480,27 @@ int main(int argc, char *argv[]) {
}
+ // output reduction factors
+ std::vector<double> avgReductionFactors(allReductionFactors.size());
+ //std::cout << "allReductionFactors.size: " << allReductionFactors.size() << std::endl;
+
+ for (size_t i=0; i<avgReductionFactors.size(); i++) {
+ avgReductionFactors[i] = 1;
+ size_t c = 0;
+
+ //std::cout << "allReductionFactors[i].size: " << allReductionFactors[i].size() << std::endl;
+ //print(allReductionFactors[i], "all reduction factors: ");
+
+ for (size_t j=0; j<allReductionFactors[i].size(); j++) {
+ avgReductionFactors[i] *= allReductionFactors[i][j];
+ c++;
+ }
+
+ avgReductionFactors[i] = std::pow(avgReductionFactors[i], 1.0/((double)c));
+ }
+
+ print(avgReductionFactors, "average reduction factors: ");
+
std::cout.rdbuf(coutbuf); //reset to standard output again
diff --git a/src/multi-body-problem/multi-body-problem.cfg b/src/multi-body-problem/multi-body-problem.cfg
index 4be02c1c..53f96e97 100644
--- a/src/multi-body-problem/multi-body-problem.cfg
+++ b/src/multi-body-problem/multi-body-problem.cfg
@@ -10,8 +10,8 @@ restarts.write = false #true
vtk.write = true
[problem]
-finalTime = 10000 # [s] #1000
-bodyCount = 2
+finalTime = 1000 # [s] #1000
+bodyCount = 3
[body]
bulkModulus = 0.5e5 # [Pa]
@@ -40,6 +40,12 @@ b = 0.017 # [ ]
a = 0.020 # [ ]
b = 0.005 # [ ]
+[boundary.neumann]
+sigmaN = 0.0 # [Pa]
+
+[boundary.dirichlet]
+finalVelocity = 5e-5 # [m/s]
+
[initialTime]
timeStep = 0
relativeTime = 0.0
@@ -47,7 +53,7 @@ relativeTau = 5e-4 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 100000
+timeSteps = 1
[u0.solver]
maximumIterations = 100
diff --git a/src/strikeslip/CMakeLists.txt b/src/strikeslip/CMakeLists.txt
new file mode 100644
index 00000000..ac71e9a9
--- /dev/null
+++ b/src/strikeslip/CMakeLists.txt
@@ -0,0 +1,43 @@
+add_custom_target(tectonic_src_strikeslip SOURCES
+ strikeslip.cfg
+ strikeslip-2D.cfg
+)
+
+set(STRIKESLIP_SOURCE_FILES
+ ../../dune/tectonic/assemblers.cc
+ ../../dune/tectonic/data-structures/body/body.cc
+ ../../dune/tectonic/data-structures/network/levelcontactnetwork.cc
+ ../../dune/tectonic/data-structures/network/contactnetwork.cc
+ ../../dune/tectonic/data-structures/enumparser.cc
+ ../../dune/tectonic/factories/strikeslipfactory.cc
+ ../../dune/tectonic/io/vtk.cc
+ ../../dune/tectonic/io/hdf5/frictionalboundary-writer.cc
+ ../../dune/tectonic/io/hdf5/iteration-writer.cc
+ #../../dune/tectonic/io/hdf5/patchinfo-writer.cc
+ ../../dune/tectonic/io/hdf5/restart-io.cc
+ ../../dune/tectonic/io/hdf5/surface-writer.cc
+ ../../dune/tectonic/io/hdf5/time-writer.cc
+ ../../dune/tectonic/problem-data/grid/simplexmanager.cc
+ ../../dune/tectonic/problem-data/grid/trianglegeometry.cc
+ ../../dune/tectonic/problem-data/grid/trianglegrids.cc
+ ../../dune/tectonic/spatial-solving/solverfactory.cc
+ ../../dune/tectonic/spatial-solving/fixedpointiterator.cc
+ ../../dune/tectonic/time-stepping/coupledtimestepper.cc
+ ../../dune/tectonic/time-stepping/adaptivetimestepper.cc
+ ../../dune/tectonic/time-stepping/rate.cc
+ ../../dune/tectonic/time-stepping/rate/rateupdater.cc
+ ../../dune/tectonic/time-stepping/state.cc
+ strikeslip.cc
+)
+
+foreach(_dim 2 3)
+ set(_strikeslip_target strikeslip-${_dim}D)
+
+ add_executable(${_strikeslip_target} ${STRIKESLIP_SOURCE_FILES})
+
+ add_dune_ug_flags(${_strikeslip_target})
+ add_dune_hdf5_flags(${_strikeslip_target})
+
+ set_property(TARGET ${_strikeslip_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+ #set_property(TARGET ${_strikeslip_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
+endforeach()
diff --git a/src/strikeslip/strikeslip-2D.cfg b/src/strikeslip/strikeslip-2D.cfg
new file mode 100644
index 00000000..05590e22
--- /dev/null
+++ b/src/strikeslip/strikeslip-2D.cfg
@@ -0,0 +1,27 @@
+# -*- mode:conf -*-
+[body0]
+smallestDiameter = 0.05 # 2e-3 [m]
+
+[body1]
+smallestDiameter = 0.05 # 2e-3 [m]
+
+[timeSteps]
+refinementTolerance = 1e-4 # 1e-5
+
+[u0.solver]
+tolerance = 1e-8
+
+[a0.solver]
+tolerance = 1e-8
+
+[v.solver]
+tolerance = 1e-8
+
+[v.fpi]
+tolerance = 1e-4 # 1e-5
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance = 1e-10
diff --git a/src/strikeslip/strikeslip.cc b/src/strikeslip/strikeslip.cc
new file mode 100644
index 00000000..412d1097
--- /dev/null
+++ b/src/strikeslip/strikeslip.cc
@@ -0,0 +1,485 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/matrices.hh>
+#include <dune/tectonic/data-structures/program_state.hh>
+#include <dune/tectonic/data-structures/friction/globalfriction.hh>
+
+#include <dune/tectonic/factories/strikeslipfactory.hh>
+
+#include <dune/tectonic/io/hdf5-writer.hh>
+#include <dune/tectonic/io/hdf5/restart-io.hh>
+#include <dune/tectonic/io/vtk.hh>
+
+#include <dune/tectonic/problem-data/bc.hh>
+#include <dune/tectonic/problem-data/mybody.hh>
+
+#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
+//#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
+#include <dune/tectonic/spatial-solving/solverfactory.hh>
+
+#include <dune/tectonic/time-stepping/adaptivetimestepper.hh>
+#include <dune/tectonic/time-stepping/rate.hh>
+#include <dune/tectonic/time-stepping/state.hh>
+#include <dune/tectonic/time-stepping/updaters.hh>
+
+#include <dune/tectonic/utils/debugutils.hh>
+#include <dune/tectonic/utils/diameter.hh>
+#include <dune/tectonic/utils/geocoordinate.hh>
+
+// for getcwd
+#include <unistd.h>
+
+//#include <tbb/tbb.h> //TODO multi threading preconditioner?
+//#include <pthread.h>
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+std::vector<std::vector<double>> allReductionFactors;
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/strikeslip/strikeslip.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/strikeslip/strikeslip-%dD.cfg", dims), parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+ return parset;
+}
+
+static std::atomic<bool> terminationRequested(false);
+void handleSignal(int signum) { terminationRequested = true; }
+
+int main(int argc, char *argv[]) {
+ try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ char buffer[256];
+ char *val = getcwd(buffer, sizeof(buffer));
+ if (val) {
+ std::cout << buffer << std::endl;
+ std::cout << argv[0] << std::endl;
+ }
+
+ std::ofstream out("strikeslip.log");
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+ auto const parset = getParameters(argc, argv);
+
+ using Assembler = MyAssembler<DefLeafGridView, dims>;
+ using field_type = Matrix::field_type;
+
+ // ----------------------
+ // set up contact network
+ // ----------------------
+ using BlocksFactory = StrikeSlipFactory<Grid, Vector>;
+ BlocksFactory blocksFactory(parset);
+ using ContactNetwork = typename BlocksFactory::ContactNetwork;
+ blocksFactory.build();
+
+ ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
+
+
+
+ const size_t bodyCount = contactNetwork.nBodies();
+
+ for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+ //printDofLocation(contactNetwork.body(i)->gridView());
+
+ //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+ //def = 1;
+ //deformedGridComplex.setDeformation(def, i);
+
+ const auto& level = *contactNetwork.level(i);
+
+ for (size_t j=0; j<level.nBodies(); j++) {
+ writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+ }
+ }
+
+ for (size_t i=0; i<bodyCount; i++) {
+ writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+ }
+
+ // ----------------------------
+ // assemble contactNetwork
+ // ----------------------------
+ contactNetwork.assemble();
+
+ //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+ // -----------------
+ // init input/output
+ // -----------------
+ std::vector<size_t> nVertices(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ nVertices[i] = contactNetwork.body(i)->nVertices();
+ }
+
+ using MyProgramState = ProgramState<Vector, ScalarVector>;
+ MyProgramState programState(nVertices);
+ auto const firstRestart = parset.get<size_t>("io.restarts.first");
+ auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
+ auto const writeRestarts = parset.get<bool>("io.restarts.write");
+ auto const writeData = parset.get<bool>("io.data.write");
+ bool const handleRestarts = writeRestarts or firstRestart > 0;
+
+
+ auto dataFile =
+ writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
+ auto restartFile = handleRestarts
+ ? std::make_unique<HDF5::File>(
+ "restarts.h5",
+ writeRestarts ? HDF5::Access::READWRITE
+ : HDF5::Access::READONLY)
+ : nullptr;
+
+
+ auto restartIO = handleRestarts
+ ? std::make_unique<RestartIO<MyProgramState>>(
+ *restartFile, nVertices)
+ : nullptr;
+
+ if (firstRestart > 0) // automatically adjusts the time and timestep
+ restartIO->read(firstRestart, programState);
+ else
+ programState.setupInitialConditions(parset, contactNetwork);
+
+ auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ for (size_t i=0; i<bodyCount; i++) {
+ contactNetwork.body(i)->setDeformation(programState.u[i]);
+ }
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();
+
+ // ------------------------
+ // assemble global friction
+ // ------------------------
+ contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+
+ auto& globalFriction = contactNetwork.globalFriction();
+ globalFriction.updateAlpha(programState.alpha);
+
+ using MyVertexBasis = typename Assembler::VertexBasis;
+ using MyCellBasis = typename Assembler::CellBasis;
+ std::vector<Vector> vertexCoordinates(bodyCount);
+ std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+ std::vector<const MyCellBasis* > cellBases(bodyCount);
+
+ auto& wPatches = blocksFactory.weakPatches();
+ std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
+
+
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ vertexBases[i] = &(body->assembler()->vertexBasis);
+ cellBases[i] = &(body->assembler()->cellBasis);
+
+ weakPatches[i].resize(1);
+ weakPatches[i][0] = wPatches[i].get();
+
+ auto& vertexCoords = vertexCoordinates[i];
+ vertexCoords.resize(nVertices[i]);
+
+ auto hostLeafView = body->grid()->hostGrid().leafGridView();
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(hostLeafView))
+ vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+ }
+
+ typename ContactNetwork::BoundaryPatches frictionBoundaries;
+ contactNetwork.boundaryPatches("friction", frictionBoundaries);
+
+ std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
+ contactNetwork.frictionPatches(frictionPatches);
+
+ auto dataWriter =
+ writeData ? std::make_unique<
+ HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+ *dataFile, vertexCoordinates, vertexBases,
+ frictionPatches, weakPatches)
+ : nullptr;
+
+ const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
+
+ IterationRegister iterationCount;
+
+ auto const report = [&](bool initial = false) {
+ if (writeData) {
+ dataWriter->reportSolution(programState, contactNetwork, globalFriction);
+ if (!initial)
+ dataWriter->reportIterations(programState, iterationCount);
+ dataFile->flush();
+ }
+
+ if (writeRestarts and !initial and
+ programState.timeStep % restartSpacing == 0) {
+ restartIO->write(programState);
+ restartFile->flush();
+ }
+
+ if (parset.get<bool>("io.printProgress"))
+ std::cout << "timeStep = " << std::setw(6) << programState.timeStep
+ << ", time = " << std::setw(12) << programState.relativeTime
+ << ", tau = " << std::setw(12) << programState.relativeTau
+ << std::endl;
+
+ if (parset.get<bool>("io.vtk.write")) {
+ std::vector<ScalarVector> stress(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(),
+ programState.u[i], stress[i]);
+
+ }
+
+ vtkWriter.write(programState.timeStep, programState.u, programState.v,
+ programState.alpha, stress);
+ }
+ };
+ report(true);
+
+ // -------------------
+ // Set up TNNMG solver
+ // -------------------
+
+ BitVector totalDirichletNodes;
+ contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+ /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+ bool val = false;
+ for (size_t d=0; d<dims; d++) {
+ val = val || totalDirichletNodes[i][d];
+ }
+
+ totalDirichletNodes[i] = val;
+ for (size_t d=0; d<dims; d++) {
+ totalDirichletNodes[i][d] = val;
+ }
+ }*/
+
+ //print(totalDirichletNodes, "totalDirichletNodes:");
+
+ //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
+ using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+ using NonlinearFactory = SolverFactory<Functional, BitVector>;
+
+ using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+ using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+ using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+ StateUpdater<ScalarVector, Vector>>;
+
+ BoundaryFunctions velocityDirichletFunctions;
+ contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+ BoundaryNodes dirichletNodes;
+ contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+ /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+ for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+ print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+ }
+ }*/
+
+ std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+ contactNetwork.frictionNodes(frictionNodes);
+
+ /*for (size_t i=0; i<frictionNodes.size(); i++) {
+ print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+ }*/
+
+ //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
+ Updaters current(
+ initRateUpdater(
+ parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunctions,
+ dirichletNodes,
+ contactNetwork.matrices(),
+ programState.u,
+ programState.v,
+ programState.a),
+ initStateUpdater<ScalarVector, Vector>(
+ parset.get<Config::stateModel>("boundary.friction.stateModel"),
+ programState.alpha,
+ nBodyAssembler.getContactCouplings(),
+ contactNetwork.couplings())
+ );
+
+
+ auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
+
+ const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+
+ auto const mustRefine = [&](Updaters &coarseUpdater,
+ Updaters &fineUpdater) {
+
+ //return false;
+ //std::cout << "Step 1" << std::endl;
+
+ std::vector<ScalarVector> coarseAlpha;
+ coarseAlpha.resize(bodyCount);
+ coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+ //print(coarseAlpha, "coarseAlpha:");
+
+ std::vector<ScalarVector> fineAlpha;
+ fineAlpha.resize(bodyCount);
+ fineUpdater.state_->extractAlpha(fineAlpha);
+
+ //print(fineAlpha, "fineAlpha:");
+
+ //std::cout << "Step 3" << std::endl;
+
+ ScalarVector::field_type energyNorm = 0;
+ for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+ //std::cout << "for " << i << std::endl;
+
+ //std::cout << not stateEnergyNorms[i] << std::endl;
+
+ if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
+ continue;
+
+ energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
+ }
+ std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
+ std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
+ return energyNorm > refinementTolerance;
+ };
+
+
+ std::signal(SIGXCPU, handleSignal);
+ std::signal(SIGINT, handleSignal);
+ std::signal(SIGTERM, handleSignal);
+
+/*
+ // set patch preconditioner for linear correction in TNNMG method
+ using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+ using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
+
+ const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+
+ auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ PatchSolver patchSolver(gsStep,
+ preconditionerParset.get<size_t>("maximumIterations"),
+ preconditionerParset.get<double>("tolerance"),
+ nullptr,
+ preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
+ false); // absolute error
+
+ Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+ Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
+ preconditioner.setPatchSolver(patchSolver);
+ preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+*/
+ // set adaptive time stepper
+ typename ContactNetwork::ExternalForces externalForces;
+ contactNetwork.externalForces(externalForces);
+
+ AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+ adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
+ programState.relativeTime, programState.relativeTau,
+ externalForces, stateEnergyNorms, mustRefine);
+
+ size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+
+ while (!adaptiveTimeStepper.reachedEnd()) {
+ programState.timeStep++;
+
+ //preconditioner.build();
+ iterationCount = adaptiveTimeStepper.advance();
+
+ programState.relativeTime = adaptiveTimeStepper.relativeTime_;
+ programState.relativeTau = adaptiveTimeStepper.relativeTau_;
+ current.rate_->extractDisplacement(programState.u);
+ current.rate_->extractVelocity(programState.v);
+ current.rate_->extractAcceleration(programState.a);
+ current.state_->extractAlpha(programState.alpha);
+ globalFriction.updateAlpha(programState.alpha);
+
+ /*print(programState.u, "current u:");
+ print(programState.v, "current v:");
+ print(programState.a, "current a:");
+ print(programState.alpha, "current alpha:");*/
+
+ contactNetwork.setDeformation(programState.u);
+
+ report();
+
+ if (programState.timeStep==timeSteps) {
+ std::cout << "limit of timeSteps reached!" << std::endl;
+ break; // TODO remove after debugging
+ }
+
+ if (terminationRequested) {
+ std::cerr << "Terminating prematurely" << std::endl;
+ break;
+ }
+
+
+ }
+
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+
+ } catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+ } catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+ }
+}
diff --git a/src/strikeslip/strikeslip.cfg b/src/strikeslip/strikeslip.cfg
new file mode 100644
index 00000000..294e3ea3
--- /dev/null
+++ b/src/strikeslip/strikeslip.cfg
@@ -0,0 +1,112 @@
+# -*- mode:conf -*-
+gravity = 9.81 # [m/s^2]
+
+[body0]
+length = 0.5 # [m]
+height = 0.5 # [m]
+depth = 0.12 # [m]
+bulkModulus = 1.5e5 # [Pa] #2190
+poissonRatio = 0.11 # [1] #0.11
+[body0.elastic]
+distFromDiag = 0.2
+density = 1300 # [kg/m^3] #750
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+[body0.viscoelastic]
+density = 1300 # [kg/m^3]
+shearViscosity = 1e4 # [Pas]
+bulkViscosity = 1e4 # [Pas]
+
+[body1]
+length = 0.5 # [m]
+height = 0.5 # [m]
+depth = 0.12 # [m]
+bulkModulus = 1.5e5 # [Pa]
+poissonRatio = 0.11 # [1]
+[body1.elastic]
+distFromDiag = 0.2
+density = 1300 # [kg/m^3]
+shearViscosity = 0 # [Pas]
+bulkViscosity = 0 # [Pas]
+[body1.viscoelastic]
+density = 1300 # [kg/m^3]
+shearViscosity = 1e4 # [Pas]
+bulkViscosity = 1e4 # [Pas]
+
+[boundary.friction]
+C = 6 # [Pa]
+mu0 = 0.48 # [ ]
+V0 = 1e-3 # [m/s]
+L = 1e-6 # [m]
+initialAlpha = 0 # [ ]
+stateModel = AgeingLaw
+frictionModel = Truncated #Regularised
+[boundary.friction.weakening]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
+[boundary.friction.strengthening]
+a = 0.054 # [ ]
+b = 0.074 # [ ]
+
+[boundary.neumann]
+sigmaN = 200.0 # [Pa]
+
+[boundary.dirichlet]
+finalVelocity = 1e-4 # [m/s]
+
+[io]
+data.write = true
+printProgress = true
+restarts.first = 0
+restarts.spacing= 20
+restarts.write = false #true
+vtk.write = true
+
+[problem]
+finalTime = 100 # [s] #1000
+bodyCount = 2
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 2e-4 # 1e-6
+
+[timeSteps]
+scheme = newmark
+timeSteps = 5000
+
+[u0.solver]
+maximumIterations = 100
+verbosity = full
+
+[a0.solver]
+maximumIterations = 100
+verbosity = full
+
+[v.solver]
+maximumIterations = 100
+verbosity = quiet
+
+[v.fpi]
+maximumIterations = 10000
+lambda = 0.5
+
+[solver.tnnmg.preconditioner]
+mode = additive
+patchDepth = 1
+maximumIterations = 2
+verbosity = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity = quiet
+
+[solver.tnnmg.main]
+pre = 1
+multi = 5 # number of multigrid steps
+post = 0
+
+
+
--
GitLab