From e106fe3c5f6fee933922b35085b994d5c1a5ffa3 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 23 Mar 2018 17:49:29 +0100
Subject: [PATCH] .
---
src/hdf5/restart-io.cc | 12 ++--
src/hdf5/restart-io.hh | 23 ++++++-
src/multi-body-problem.cc | 2 +-
src/program_state.hh | 133 +++++++++++++++++++++++++-------------
4 files changed, 117 insertions(+), 53 deletions(-)
diff --git a/src/hdf5/restart-io.cc b/src/hdf5/restart-io.cc
index c5724e55..0d16b86c 100644
--- a/src/hdf5/restart-io.cc
+++ b/src/hdf5/restart-io.cc
@@ -5,15 +5,15 @@
#include "restart-io.hh"
template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
- : displacementWriter_(group, "displacement", vertexCount,
+RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, const std::vector<size_t>& vertexCounts)
+ : displacementWriter_(group, "displacement", vertexCounts,
Vector::block_type::dimension),
- velocityWriter_(group, "velocity", vertexCount,
+ velocityWriter_(group, "velocity", vertexCounts,
Vector::block_type::dimension),
- accelerationWriter_(group, "acceleration", vertexCount,
+ accelerationWriter_(group, "acceleration", vertexCounts,
Vector::block_type::dimension),
- stateWriter_(group, "state", vertexCount),
- weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
+ stateWriter_(group, "state", vertexCounts),
+ weightedNormalStressWriter_(group, "weightedNormalStress", vertexCounts),
relativeTimeWriter_(group, "relativeTime"),
relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
diff --git a/src/hdf5/restart-io.hh b/src/hdf5/restart-io.hh
index f90b01a7..1cf5167a 100644
--- a/src/hdf5/restart-io.hh
+++ b/src/hdf5/restart-io.hh
@@ -4,12 +4,33 @@
#include <dune/fufem/hdf5/file.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
+template <class ProgramState> class RestartBodyIO {
+ using ScalarVector = typename ProgramState::ScalarVector;
+ using Vector = typename ProgramState::Vector;
+
+public:
+ RestartBodyIO(HDF5::Grouplike &group, const std::vector<size_t>& vertexCounts);
+
+ void write(ProgramState const &programState);
+
+ void read(size_t timeStep, ProgramState &programState);
+
+private:
+ HDF5::SequenceIO<2> displacementWriter_;
+ HDF5::SequenceIO<2> velocityWriter_;
+ HDF5::SequenceIO<2> accelerationWriter_;
+ HDF5::SequenceIO<1> stateWriter_;
+ HDF5::SequenceIO<1> weightedNormalStressWriter_;
+ HDF5::SequenceIO<0> relativeTimeWriter_;
+ HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
+};
+
template <class ProgramState> class RestartIO {
using ScalarVector = typename ProgramState::ScalarVector;
using Vector = typename ProgramState::Vector;
public:
- RestartIO(HDF5::Grouplike &group, size_t vertexCount);
+ RestartIO(HDF5::Grouplike &group, const std::vector<size_t>& vertexCounts);
void write(ProgramState const &programState);
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 54e76531..38971c5b 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -211,7 +211,7 @@ int main(int argc, char *argv[]) {
std::vector<std::shared_ptr<DefLeafGridView>> leafViews(deformedGrids.size());
std::vector<std::shared_ptr<DefLevelGridView>> levelViews(deformedGrids.size());
- std::vector<int> leafVertexCounts(deformedGrids.size());
+ std::vector<size_t> leafVertexCounts(deformedGrids.size());
using LeafFaces = MyFaces<DefLeafGridView>;
using LevelFaces = MyFaces<DefLevelGridView>;
diff --git a/src/program_state.hh b/src/program_state.hh
index 73f9bc2e..b48103cf 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -17,29 +17,62 @@
#include "matrices.hh"
#include "spatial-solving/solverfactory.hh"
-template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
+template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
public:
using Vector = VectorTEMPLATE;
using ScalarVector = ScalarVectorTEMPLATE;
- ProgramState(const std::vector<int>& leafVertexCounts)
- : bodyCount(leafVertexCounts.size()) {
- u.resize(bodyCount);
- v.resize(bodyCount);
- a.resize(bodyCount);
- alpha.resize(bodyCount);
- weightedNormalStress.resize(bodyCount);
-
- for (size_t i=0; i<bodyCount; i++) {
- size_t leafVertexCount = leafVertexCounts[i];
- u[i].resize(leafVertexCount);
- v[i].resize(leafVertexCount);
- a[i].resize(leafVertexCount);
- alpha[i].resize(leafVertexCount);
- weightedNormalStress[i].resize(leafVertexCount);
+ BodyState() {}
+
+ BodyState(Vector * _u, Vector * _v, Vector * _a, ScalarVector * _alpha, ScalarVector * _weightedNormalStress)
+ : u(_u),
+ v(_v),
+ a(_a),
+ alpha(_alpha),
+ weightedNormalStress(_weightedNormalStress) {}
+
+public:
+ Vector * const u;
+ Vector * const v;
+ Vector * const a;
+ ScalarVector * const alpha;
+ ScalarVector * const weightedNormalStress;
+};
+
+
+template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
+public:
+ using Vector = VectorTEMPLATE;
+ using ScalarVector = ScalarVectorTEMPLATE;
+ using BodyState = BodyState<Vector, ScalarVector>;
+
+ ProgramState(const std::vector<size_t>& leafVertexCounts)
+ : bodyCount_(leafVertexCounts.size()),
+ bodies(bodyCount_),
+ u(bodyCount_),
+ v(bodyCount_),
+ a(bodyCount_),
+ alpha(bodyCount_),
+ weightedNormalStress(bodyCount_) {
+
+ for (size_t i=0; i<bodyCount_; i++) {
+ size_t leafVertexCount = leafVertexCounts[i];
+
+ u[i].resize(leafVertexCount),
+ v[i].resize(leafVertexCount),
+ a[i].resize(leafVertexCount),
+ alpha[i].resize(leafVertexCount),
+ weightedNormalStress[i].resize(leafVertexCount),
+
+ bodies[i] = BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
}
}
+ size_t size() {
+ return bodyCount_;
+ }
+
+
// Set up initial conditions
template <class Matrix, class GridView>
void setupInitialConditions(
@@ -57,32 +90,48 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
using LocalMatrix = typename Matrix::block_type;
auto constexpr dims = LocalVector::dimension;
- /*
+
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
- Dune::BitSetVector<dims> const &_dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrix,
- Vector const &_rhs, Vector &_x,
+ const std::vector<Dune::BitSetVector<dims>>& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+ const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
Dune::ParameterTree const &_localParset) {
- using LinearFactory = SolverFactory<
- dims, BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- typename GridView::Grid>;
- ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
- LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
- assemblers.gridView.grid(), _dirichletNodes);
+ std::vector<const Matrix*> matrices_ptr(_matrices.size());
+ for (size_t i=0; i<matrices_ptr.size(); i++) {
+ matrices_ptr[i] = _matrices[i].get();
+ }
- typename LinearFactory::ConvexProblem convexProblem(
- 1.0, _matrix, zeroNonlinearity, _rhs, _x);
- typename LinearFactory::BlockProblem problem(parset, convexProblem);
+ /*std::vector<Matrix> matrices(velocityMatrices.size());
+ std::vector<Vector> rhs(velocityRHSs.size());
+ for (size_t i=0; i<globalFriction_.size(); i++) {
+ matrices[i] = velocityMatrices[i];
+ rhs[i] = velocityRHSs[i];
- auto multigridStep = factory.getStep();
- multigridStep->setProblem(_x, problem);
- //multigridStep->setProblem(_x);
- EnergyNorm<Matrix, Vector> const norm(_matrix);
+ globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
+ globalFriction_[i]->addGradient(v_rel[i], rhs[i]);
+
+ matrices_ptr[i] = &matrices[i];
+ }*/
+
+ // assemble full global contact problem
+ Matrix bilinearForm;
+ nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs;
+ nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
+ Vector totalX;
+ nBodyAssembler.nodalToTransformed(_x, totalX);
+ using LinearFactory = SolverFactory<typename GridView::Grid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+ LinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
+
+ auto multigridStep = factory.getStep();
+ multigridStep->setProblem(bilinearForm, totalX, totalRhs);
+
+ const EnergyNorm<Matrix, Vector> norm(bilinearForm);
LoopSolver<Vector> solver(
multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
@@ -93,21 +142,15 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
solver.preprocess();
solver.solve();
- Vector totalX = multigridStep->getSol();
-
- // cleanup
- delete(multigridStep);
-
- nBodyAssembler.postprocess(totalX, x);
+ nBodyAssembler.postprocess(multigridStep->getSol(), _x);
};
- */
timeStep = 0;
relativeTime = 0.0;
relativeTau = 1e-6;
- std::vector<Vector> ell0(bodyCount);
- for (size_t i=0; i<bodyCount; i++) {
+ std::vector<Vector> ell0(bodyCount_);
+ for (size_t i=0; i<bodyCount_; i++) {
// Initial velocity
v[i] = 0.0;
@@ -124,7 +167,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Initial acceleration: Computed in agreement with Ma = ell0 - Au
// (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
std::vector<Vector> accelerationRHS = ell0;
- for (size_t i=0; i<bodyCount; i++) {
+ for (size_t i=0; i<bodyCount_; i++) {
// Initial state
alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
@@ -141,6 +184,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
}
public:
+ std::vector<const BodyState> bodies;
std::vector<Vector> u;
std::vector<Vector> v;
std::vector<Vector> a;
@@ -151,7 +195,6 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
size_t timeStep;
private:
- const size_t bodyCount;
-};
+ const size_t bodyCount_;
#endif
--
GitLab