Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • podlesny/dune-tectonic
  • agnumpde/dune-tectonic
2 results
Show changes
Showing
with 757 additions and 94 deletions
#ifndef SRC_CONTACTNETWORKFACTORY_HH
#define SRC_CONTACTNETWORKFACTORY_HH
#include <dune/common/parametertree.hh>
#include "../data-structures/contactnetwork.hh"
#include "../data-structures/levelcontactnetwork.hh"
#include "levelcontactnetworkfactory.hh"
template <class GridType, int dim>
class ContactNetworkFactory {
public:
using ContactNetwork = ContactNetwork<GridType, dim>;
protected:
using LevelContactNetworkFactory = LevelContactNetworkFactory<GridType, dim>;
const Dune::ParameterTree& parset_;
std::vector<std::shared_ptr<LevelContactNetworkFactory>> levelContactNetworkFactories_;
ContactNetwork contactNetwork_;
private:
virtual void setBodies() = 0;
virtual void setCouplings() = 0;
virtual void setBoundaryConditions() = 0;
public:
ContactNetworkFactory(const Dune::ParameterTree& parset) :
parset_(parset),
contactNetwork_() {}
void build() {
setBodies();
contactNetwork_.setBodies();
setCouplings();
contactNetwork_.setCouplings();
setBoundaryConditions();
}
ContactNetwork& levelContactNetwork() {
return contactNetwork_;
}
};
#endif
#ifndef SRC_LEVELCONTACTNETWORKFACTORY_HH
#define SRC_LEVELCONTACTNETWORKFACTORY_HH
#include <dune/common/parametertree.hh>
#include "../data-structures/levelcontactnetwork.hh"
template <class GridType, int dims>
class LevelContactNetworkFactory {
public:
using LevelContactNetwork = LevelContactNetwork<GridType, dims>;
protected:
using Body = typename LevelContactNetwork::Body;
using FrictionCouplingPair = typename LevelContactNetwork::FrictionCouplingPair;
const Dune::ParameterTree& parset_;
const size_t bodyCount_;
const size_t couplingCount_;
std::vector<std::shared_ptr<Body>> bodies_;
std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
LevelContactNetwork levelContactNetwork_;
private:
virtual void setBodies() = 0;
virtual void setCouplings() = 0;
virtual void setBoundaryConditions() = 0;
public:
LevelContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
parset_(parset),
bodyCount_(bodyCount),
couplingCount_(couplingCount),
bodies_(bodyCount_),
couplings_(couplingCount_),
levelContactNetwork_(bodyCount_, couplingCount_) {}
void build() {
setBodies();
levelContactNetwork_.setBodies(bodies_);
setCouplings();
levelContactNetwork_.setCouplings(couplings_);
setBoundaryConditions();
}
LevelContactNetwork& levelContactNetwork() {
return levelContactNetwork_;
}
};
#endif
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/contact/projections/normalprojection.hh>
#include "../multi-body-problem-data/bc.hh"
#include "../multi-body-problem-data/grid/cuboidgeometry.hh"
#include "../multi-body-problem-data/myglobalfrictiondata.hh"
#include "../multi-body-problem-data/weakpatch.hh"
#include "../utils/diameter.hh"
#include "stackedblocksfactory.hh"
template <class GridType, int dims>
void StackedBlocksFactory<GridType, dims>::setBodies() {
// set up cuboid geometries
double const length = 1.00;
double const width = 0.27;
double const weakLength = 0.20;
std::vector<LocalVector> origins(this->bodyCount_);
std::vector<LocalVector> lowerWeakOrigins(this->bodyCount_);
std::vector<LocalVector> upperWeakOrigins(this->bodyCount_);
#if MY_DIM == 3
double const depth = 0.60;
origins[0] = {0,0,0};
lowerWeakOrigins[0] = {0.2,0,0};
upperWeakOrigins[0] = {0.2,width,0};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width, depth);
for (size_t i=1; i<this->bodyCount_-1; i++) {
origins[i] = cuboidGeometries_[i-1]->D;
lowerWeakOrigins[i] = {0.6,i*width,0};
upperWeakOrigins[i] = {0.6,(i+1)*width,0};
cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width, depth);
}
const size_t idx = this->bodyCount_-1;
origins[idx] = cuboidGeometries_[idx-1]->D;
lowerWeakOrigins[idx] = {0.6,idx*width,0};
upperWeakOrigins[idx] = {0.6,(idx+1)*width,0};
cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width, depth);
#elif MY_DIM == 2
origins[0] = {0,0};
lowerWeakOrigins[0] = {0.2,0};
upperWeakOrigins[0] = {0.6,width};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width);
for (size_t i=1; i<this->bodyCount_-1; i++) {
origins[i] = cuboidGeometries_[i-1]->D;
lowerWeakOrigins[i] = {0.6,i*width};
upperWeakOrigins[i] = {0.6,(i+1)*width};
cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width);
}
const size_t idx = this->bodyCount_-1;
origins[idx] = cuboidGeometries_[idx-1]->D + LocalVector({0.25,0});
lowerWeakOrigins[idx] = {0.6,idx*width};
upperWeakOrigins[idx] = {0.6,(idx+1)*width};
cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, 5.0/8*length, width);
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
// set up reference grids
gridConstructor_ = new GridsConstructor<GridType>(cuboidGeometries_);
auto& grids = gridConstructor_->getGrids();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
// define weak patch and refine grid
auto lowerWeakPatch = lowerWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
auto upperWeakPatch = upperWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
getWeakPatch<LocalVector>(this->parset_.sub("boundary.friction.weakening"), cuboidGeometry, *lowerWeakPatch, *upperWeakPatch);
refine(*grids[i], *lowerWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
refine(*grids[i], *upperWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
// determine minDiameter and maxDiameter
double minDiameter = std::numeric_limits<double>::infinity();
double maxDiameter = 0.0;
for (auto &&e : elements(grids[i]->leafGridView())) {
auto const geometry = e.geometry();
auto const diam = diameter(geometry);
minDiameter = std::min(minDiameter, diam);
maxDiameter = std::max(maxDiameter, diam);
}
std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
}
for (size_t i=0; i<this->bodyCount_; i++) {
this->bodies_[i] = std::make_shared<typename Base::Body>(bodyData_, grids[i]);
}
}
template <class GridType, int dims>
void StackedBlocksFactory<GridType, dims>::setCouplings() {
const auto deformedGrids = this->levelContactNetwork_.deformedGrids();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
leafFaces_[i] = std::make_shared<LeafFaces>(this->levelContactNetwork_.leafView(i), cuboidGeometry);
levelFaces_[i] = std::make_shared<LevelFaces>(this->levelContactNetwork_.levelView(i), cuboidGeometry);
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
for (size_t i=0; i<this->couplings_.size(); i++) {
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);
coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
coupling->setWeakeningPatch(upperWeakPatches_[i]);
coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], CuboidGeometry::lengthScale));
}
}
template <class GridType, int dims>
void StackedBlocksFactory<GridType, dims>::setBoundaryConditions() {
using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();
const double lengthScale = CuboidGeometry::lengthScale;
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& body = this->levelContactNetwork_.body(i);
const auto& leafVertexCount = body->leafVertexCount();
std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
// Neumann boundary
std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
body->addBoundaryCondition(neumannBoundary);
// right Dirichlet Boundary
// identify Dirichlet nodes on leaf level
std::shared_ptr<Dune::BitSetVector<dims>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(leafVertexCount);
for (int j=0; j<leafVertexCount; j++) {
if (leafFaces_[i]->upper.containsVertex(j))
(*velocityDirichletNodes)[j][0] = true;
#if MY_DIM == 3 //TODO: wrong, needs revision
if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
zeroDirichletNodes->at(j)[2] = true;
#endif
}
if (i>0) {
std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
velocityDirichletBoundary->setBoundaryPatch(body->leafView(), velocityDirichletNodes);
velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
body->addBoundaryCondition(velocityDirichletBoundary);
} else {
//velocityDirichletBoundary->setBoundaryFunction(neumannFunction);
}
}
// lower Dirichlet Boundary
const auto& firstBody = this->levelContactNetwork_.body(0);
const auto& firstLeafVertexCount = firstBody->leafVertexCount();
std::shared_ptr<Dune::BitSetVector<dims>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dims>>(firstLeafVertexCount);
for (int j=0; j<firstLeafVertexCount; j++) {
if (leafFaces_[0]->lower.containsVertex(j)) {
for (size_t d=0; d<dims; d++) {
(*zeroDirichletNodes)[j][d] = true;
}
}
}
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
zeroDirichletBoundary->setBoundaryPatch(firstBody->leafView(), zeroDirichletNodes);
std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
firstBody->addBoundaryCondition(zeroDirichletBoundary);
}
#include "stackedblocksfactory_tmpl.cc"
#ifndef SRC_STACKEDBLOCKSFACTORY_HH
#define SRC_STACKEDBLOCKSFACTORY_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/fufem/boundarypatch.hh>
#include "levelcontactnetworkfactory.hh"
#include "../multi-body-problem-data/mybody.hh"
#include "../multi-body-problem-data/grid/mygrids.hh"
template <class GridType, int dims> class StackedBlocksFactory : public LevelContactNetworkFactory<GridType, dims>{
private:
using Base = LevelContactNetworkFactory<GridType, dims>;
using LocalVector = typename Base::LevelContactNetwork::LocalVector;
using DeformedLeafGridView = typename Base::LevelContactNetwork::DeformedLeafGridView;
using DeformedLevelGridView = typename Base::LevelContactNetwork::DeformedLevelGridView;
using LevelBoundaryPatch = BoundaryPatch<DeformedLevelGridView>;
using LeafBoundaryPatch = BoundaryPatch<DeformedLeafGridView>;
using LeafFaces = MyFaces<DeformedLeafGridView>;
using LevelFaces = MyFaces<DeformedLevelGridView>;
private:
const std::shared_ptr<MyBodyData<dims>> bodyData_; // material properties of bodies
GridsConstructor<GridType>* gridConstructor_;
std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> lowerWeakPatches_;
std::vector<std::shared_ptr<ConvexPolyhedron<LocalVector>>> upperWeakPatches_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
public:
StackedBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
bodyData_(std::make_shared<MyBodyData<dims>>(this->parset_, zenith_())),
cuboidGeometries_(this->bodyCount_),
leafFaces_(this->bodyCount_),
levelFaces_(this->bodyCount_),
lowerWeakPatches_(this->bodyCount_),
upperWeakPatches_(this->bodyCount_),
nonmortarPatch_(this->couplingCount_),
mortarPatch_(this->couplingCount_)
{}
~StackedBlocksFactory() {
delete gridConstructor_;
}
void setBodies();
void setCouplings();
void setBoundaryConditions();
private:
static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
#if MY_DIM == 2
return {0, 1};
#elif MY_DIM == 3
return {0, 1, 0};
#endif
}
};
#endif
#ifndef MY_DIM
#error MY_DIM unset
#endif
#include "../explicitgrid.hh"
template class StackedBlocksFactory<Grid, MY_DIM>;
#ifndef SRC_FRICTIONCOUPLINGPAIR_HH
#define SRC_FRICTIONCOUPLINGPAIR_HH
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/contact/common/couplingpair.hh>
#include <dune/tectonic/globalfrictiondata.hh>
template <class GridType, class LocalVectorType, class field_type = double>
class FrictionCouplingPair : public Dune::Contact::CouplingPair<GridType, GridType, field_type>{
private:
static const int dim = GridType::dimensionworld;
using Base = Dune::Contact::CouplingPair<GridType,GridType,field_type>;
using LocalVector = LocalVectorType;
// friction data
std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
std::shared_ptr<GlobalFrictionData<dim>> frictionData_;
public:
void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
weakeningPatch_ = weakeningPatch;
}
void setFrictionData(std::shared_ptr<GlobalFrictionData<dim>> frictionData) {
frictionData_ = frictionData;
}
const auto& weakeningPatch() const {
return *weakeningPatch_;
}
const auto& frictionData() const {
return *frictionData_;
}
};
#endif
#ifndef SRC_HDF5_WRITER_HH
#define SRC_HDF5_WRITER_HH
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/hdf5/file.hh>
#include "hdf5/frictionalboundary-writer.hh"
#include "hdf5/iteration-writer.hh"
#include "hdf5/patchinfo-writer.hh"
#include "hdf5/surface-writer.hh"
#include "hdf5/time-writer.hh"
template <class ProgramState, class VertexBasis, class GridView>
class HDF5Writer {
private:
using Vector = typename ProgramState::Vector;
using Patch = BoundaryPatch<GridView>;
using LocalVector = typename Vector::block_type;
public:
HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates,
VertexBasis const &vertexBasis, Patch const &surface,
Patch const &frictionalBoundary,
ConvexPolyhedron<LocalVector> const &weakPatch)
: file_(file),
iterationWriter_(file_),
timeWriter_(file_),
#if MY_DIM == 3
patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
#endif
surfaceWriter_(file_, vertexCoordinates, surface),
frictionalBoundaryWriter_(file_, vertexCoordinates,
frictionalBoundary) {
}
template <class Friction>
void reportSolution(ProgramState const &programState,
// for the friction coefficient
Friction &friction) {
timeWriter_.write(programState);
#if MY_DIM == 3
patchInfoWriter_.write(programState);
#endif
surfaceWriter_.write(programState);
frictionalBoundaryWriter_.write(programState, friction);
}
void reportIterations(ProgramState const &programState,
IterationRegister const &iterationCount) {
iterationWriter_.write(programState.timeStep, iterationCount);
}
private:
HDF5::Grouplike &file_;
IterationWriter iterationWriter_;
TimeWriter<ProgramState> timeWriter_;
#if MY_DIM == 3
PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
#endif
SurfaceWriter<ProgramState, GridView> surfaceWriter_;
FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
};
#endif
#include "../explicitvectors.hh"
#include "../explicitgrid.hh"
#include "../program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>;
using P1Basis = P1NodalBasis<GridView, double>;
using MyFunction = BasisGridFunction<P1Basis, Vector>;
template class GridEvaluator<LocalVector, GridView>;
template Dune::Matrix<typename MyFunction::RangeType>
GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;
#ifndef SRC_HDF5_BODYWRITER_HH
#define SRC_HDF5_BODYWRITER_HH
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/hdf5/file.hh>
#include "hdf5/frictionalboundary-writer.hh"
#include "hdf5/patchinfo-writer.hh"
template <class ProgramState, class VertexBasis, class GridView>
class HDF5BodyWriter {
private:
using VertexCoordinates = typename ProgramState::Vector;
using FrictionPatches = std::vector<const BoundaryPatch<GridView>* >;
using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
using LocalVector = typename VertexCoordinates::block_type;
friend class HDF5LevelWriter<ProgramState, VertexBasis, GridView>;
public:
HDF5BodyWriter(
HDF5::Grouplike &file,
const VertexCoordinates& vertexCoordinates,
const VertexBasis& vertexBasis,
const FrictionPatches& frictionPatches,
const WeakPatches& weakPatches) :
patchCount_(frictionPatches.size()),
file_(file),
#if MY_DIM == 3
patchInfoWriters_(patchCount_),
#endif
frictionBoundaryWriters_(patchCount_) {
assert(patchCount_ == weakPatches.size());
for (size_t i=0; i<patchCount_; i++) {
#if MY_DIM == 3
patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, vertexBasis, *frictionPatches[i], *weakPatches[i], i);
#endif
frictionBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates, *frictionPatches[i], i);
}
}
~HDF5BodyWriter() {
for (size_t i=0; i<patchCount_; i++) {
#if MY_DIM == 3
delete patchInfoWriters_[i];
#endif
delete frictionBoundaryWriters_[i];
}
}
template <class Friction>
void reportSolution(ProgramState const &programState,
// for the friction coefficient
std::vector<std::shared_ptr<Friction>>& friction) {
for (size_t i=0; i<patchCount_; i++) {
#if MY_DIM == 3
patchInfoWriters_[i]->write(programState);
#endif
frictionBoundaryWriters_[i]->write(programState, *friction[i]);
}
}
private:
const size_t patchCount_;
HDF5::Grouplike &file_;
#if MY_DIM == 3
std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_;
#endif
std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionBoundaryWriters_;
};
#endif
#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 {
private:
using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
using VertexBases = std::vector<VertexBasis>;
using FrictionPatches = std::vector<typename HDF5BodyWriter::FrictionPatches>;
using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
public:
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] = new HDF5BodyWriter<ProgramState, VertexBasis, GridView>(file_, vertexCoordinates[i], vertexBases[i], frictionBoundaries[i], weakPatches[i]);
}
}
~HDF5LevelWriter() {
for (size_t i=0; i<bodyCount_; i++) {
delete bodyWriters_[i];
}
}
template <class Friction>
void reportSolution(
ProgramState const &programState,
// for the friction coefficient
std::vector<std::shared_ptr<Friction>>& friction) {
timeWriter_.write(programState);
for (size_t i=0; i<bodyCount_; i++) {
bodyWriters_[i]->reportSolution(programState, *friction[i]); //TODO
}
}
void reportIterations(
ProgramState const &programState,
IterationRegister const &iterationCount) {
iterationWriter_.write(programState.timeStep, iterationCount);
}
private:
const size_t bodyCount_;
HDF5::Grouplike &file_;
IterationWriter iterationWriter_;
TimeWriter<ProgramState> timeWriter_;
std::vector<HDF5BodyWriter* > bodyWriters_;
};
#endif
......@@ -10,8 +10,9 @@
template <class ProgramState, class GridView>
FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
HDF5::Grouplike &file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary)
: group_(file, "frictionalBoundary"),
Patch const &frictionalBoundary, size_t const id)
: id_(id),
group_(file, "frictionalBoundary"+std::to_string(id_)),
frictionalBoundary_(frictionalBoundary),
frictionalBoundaryDisplacementWriter_(group_, "displacement",
frictionalBoundary.numVertices(),
......@@ -35,23 +36,24 @@ template <class ProgramState, class GridView>
template <class Friction>
void FrictionalBoundaryWriter<ProgramState, GridView>::write(
ProgramState const &programState, Friction &friction) {
auto const frictionalBoundaryDisplacements =
restrictToSurface(programState.u, frictionalBoundary_);
restrictToSurface(programState.u[id_], frictionalBoundary_);
addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
frictionalBoundaryDisplacements);
auto const frictionalBoundaryVelocities =
restrictToSurface(programState.v, frictionalBoundary_);
restrictToSurface(programState.v[id_], frictionalBoundary_);
addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
frictionalBoundaryVelocities);
auto const frictionalBoundaryStates =
restrictToSurface(programState.alpha, frictionalBoundary_);
restrictToSurface(programState.alpha[id_], frictionalBoundary_);
addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
frictionalBoundaryStates);
friction.updateAlpha(programState.alpha);
auto const c = friction.coefficientOfFriction(programState.v);
auto const c = friction.coefficientOfFriction(programState.v[id_]);
auto const frictionalBoundaryCoefficient =
restrictToSurface(c, frictionalBoundary_);
addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
......
......@@ -11,13 +11,15 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
using Patch = BoundaryPatch<GridView>;
public:
FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary);
FrictionalBoundaryWriter(HDF5::Grouplike& file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary, size_t const id);
template <class Friction>
void write(ProgramState const &programState, Friction &friction);
private:
size_t const id_;
HDF5::Group group_;
Patch const &frictionalBoundary_;
......
#include "../explicitvectors.hh"
#include "../explicitgrid.hh"
#include "../../explicitvectors.hh"
#include "../../explicitgrid.hh"
#include "../program_state.hh"
#include "../../data-structures/program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>;
using MyFriction = GlobalFriction<Matrix, Vector>;
template class FrictionalBoundaryWriter<MyProgramState, GridView>;
template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
template class FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>;
template void FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>::write(
MyProgramState const &programState, MyFriction &friction);
......@@ -3,7 +3,7 @@
#include <dune/fufem/hdf5/sequenceio.hh>
#include "../time-stepping/adaptivetimestepper.hh"
#include "../../time-stepping/adaptivetimestepper.hh"
class IterationWriter {
public:
......
......@@ -68,8 +68,10 @@ template <class ProgramState, class VertexBasis, class GridView>
PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
HDF5::Grouplike &file, VertexBasis const &vertexBasis,
Patch const &frictionalBoundary,
ConvexPolyhedron<LocalVector> const &weakPatch)
: group_(file, "weakPatchGrid"),
ConvexPolyhedron<LocalVector> const &weakPatch,
size_t const id)
: id_(id),
group_(file, "weakPatchGrid"+std::to_string(id_)),
vertexBasis_(vertexBasis),
gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
weakPatchGridVelocityWriter_(
......@@ -87,7 +89,8 @@ PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
template <class ProgramState, class VertexBasis, class GridView>
void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
ProgramState const &programState) {
BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
// TODO
BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]);
auto const gridVelocity = gridEvaluator_.evaluate(velocity);
addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
}
......
......@@ -8,7 +8,7 @@
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
#include "../one-body-problem-data/mygeometry.hh"
#include "../../one-body-problem-data/mygeometry.hh"
template <class LocalVector, class GridView> class GridEvaluator {
using Element = typename GridView::Grid::template Codim<0>::Entity;
......@@ -36,11 +36,14 @@ class PatchInfoWriter {
public:
PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
Patch const &frictionalBoundary,
ConvexPolyhedron<LocalVector> const &weakPatch);
ConvexPolyhedron<LocalVector> const &weakPatch,
size_t const id);
void write(ProgramState const &programState);
private:
const size_t id_;
HDF5::Group group_;
VertexBasis const &vertexBasis_;
......
#include "../../explicitvectors.hh"
#include "../../explicitgrid.hh"
#include "../../data-structures/program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>;
using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
using MyFunction = BasisGridFunction<P1Basis, Vector>;
template class GridEvaluator<LocalVector, DeformedGrid::LeafGridView>;
template Dune::Matrix<typename MyFunction::RangeType>
GridEvaluator<LocalVector, DeformedGrid::LeafGridView>::evaluate(MyFunction const &f) const;
template class PatchInfoWriter<MyProgramState, P1Basis, DeformedGrid::LeafGridView>;
......@@ -5,26 +5,66 @@
#include "restart-io.hh"
template <class ProgramState>
RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
: displacementWriter_(group, "displacement", vertexCount,
Vector::block_type::dimension),
velocityWriter_(group, "velocity", vertexCount,
Vector::block_type::dimension),
accelerationWriter_(group, "acceleration", vertexCount,
Vector::block_type::dimension),
stateWriter_(group, "state", vertexCount),
weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
RestartBodyIO<ProgramState>::RestartBodyIO(HDF5::Grouplike &group, const size_t vertexCount, const size_t id)
: id_(id),
displacementWriter_(group, "displacement"+std::to_string(id_), vertexCount,
ProgramState::Vector::block_type::dimension),
velocityWriter_(group, "velocity"+std::to_string(id_), vertexCount,
ProgramState::Vector::block_type::dimension),
accelerationWriter_(group, "acceleration"+std::to_string(id_), vertexCount,
ProgramState::Vector::block_type::dimension),
stateWriter_(group, "state"+std::to_string(id_), vertexCount),
weightedNormalStressWriter_(group, "weightedNormalStress"+std::to_string(id_), vertexCount) {}
template <class ProgramState>
void RestartBodyIO<ProgramState>::write(const ProgramState& programState) {
addEntry(displacementWriter_, programState.timeStep, programState.u[id_]);
addEntry(velocityWriter_, programState.timeStep, programState.v[id_]);
addEntry(accelerationWriter_, programState.timeStep, programState.a[id_]);
addEntry(stateWriter_, programState.timeStep, programState.alpha[id_]);
addEntry(weightedNormalStressWriter_, programState.timeStep,
programState.weightedNormalStress[id_]);
}
template <class ProgramState>
void RestartBodyIO<ProgramState>::read(size_t timeStep,
ProgramState& programState) {
readEntry(displacementWriter_, timeStep, programState.u[id_]);
readEntry(velocityWriter_, timeStep, programState.v[id_]);
readEntry(accelerationWriter_, timeStep, programState.a[id_]);
readEntry(stateWriter_, timeStep, programState.alpha[id_]);
readEntry(weightedNormalStressWriter_, timeStep,
programState.weightedNormalStress[id_]);
}
template <class ProgramState>
RestartIO<ProgramState>::RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts)
: bodyIO_(vertexCounts.size()),
relativeTimeWriter_(group, "relativeTime"),
relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {
for (size_t i=0; i<bodyIO_.size(); i++) {
bodyIO_[i] = new RestartBodyIO<ProgramState>(group, vertexCounts[i], i);
}
}
template <class ProgramState>
RestartIO<ProgramState>::~RestartIO() {
for (size_t i=0; i<bodyIO_.size(); i++) {
delete bodyIO_[i];
}
}
template <class ProgramState>
void RestartIO<ProgramState>::write(ProgramState const &programState) {
addEntry(displacementWriter_, programState.timeStep, programState.u);
addEntry(velocityWriter_, programState.timeStep, programState.v);
addEntry(accelerationWriter_, programState.timeStep, programState.a);
addEntry(stateWriter_, programState.timeStep, programState.alpha);
addEntry(weightedNormalStressWriter_, programState.timeStep,
programState.weightedNormalStress);
assert(programState.size() == bodyIO_.size());
for (size_t i=0; i<bodyIO_.size(); i++) {
bodyIO_[i]->write(programState);
}
addEntry(relativeTimeWriter_, programState.timeStep,
programState.relativeTime);
addEntry(relativeTimeIncrementWriter_, programState.timeStep,
......@@ -33,14 +73,15 @@ void RestartIO<ProgramState>::write(ProgramState const &programState) {
template <class ProgramState>
void RestartIO<ProgramState>::read(size_t timeStep,
ProgramState &programState) {
ProgramState& programState) {
assert(programState.size() == bodyIO_.size());
programState.timeStep = timeStep;
readEntry(displacementWriter_, timeStep, programState.u);
readEntry(velocityWriter_, timeStep, programState.v);
readEntry(accelerationWriter_, timeStep, programState.a);
readEntry(stateWriter_, timeStep, programState.alpha);
readEntry(weightedNormalStressWriter_, timeStep,
programState.weightedNormalStress);
for (size_t i=0; i<bodyIO_.size(); i++) {
bodyIO_[i]->read(timeStep, programState);
}
readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
}
......
......@@ -4,23 +4,37 @@
#include <dune/fufem/hdf5/file.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
template <class ProgramState> class RestartIO {
using ScalarVector = typename ProgramState::ScalarVector;
using Vector = typename ProgramState::Vector;
template <class ProgramState> class RestartBodyIO {
public:
RestartIO(HDF5::Grouplike &group, size_t vertexCount);
RestartBodyIO(HDF5::Grouplike& group, const size_t vertexCount, const size_t id);
void write(ProgramState const &programState);
void write(const ProgramState& programState);
void read(size_t timeStep, ProgramState &programState);
void read(size_t timeStep, ProgramState& programState);
private:
const size_t id_;
HDF5::SequenceIO<2> displacementWriter_;
HDF5::SequenceIO<2> velocityWriter_;
HDF5::SequenceIO<2> accelerationWriter_;
HDF5::SequenceIO<1> stateWriter_;
HDF5::SequenceIO<1> weightedNormalStressWriter_;
};
template <class ProgramState> class RestartIO {
public:
RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts);
~RestartIO();
void write(const ProgramState& programState);
void read(size_t timeStep, ProgramState& programState);
private:
std::vector<RestartBodyIO<ProgramState>* > bodyIO_;
HDF5::SequenceIO<0> relativeTimeWriter_;
HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
};
......