Skip to content
Snippets Groups Projects
Commit b660cb8d authored by podlesny's avatar podlesny
Browse files

write data and restarts

parent 8f824a8a
Branches
No related tags found
No related merge requests found
Showing
with 269 additions and 373 deletions
import ConfigParser as cp
import os
import numpy as np
import csv
import h5py
from support.find_quakes import find_quakes
NBODIES = 2
FINAL_TIME = 1000 # s
FINAL_VELOCITY = 1e-5 # m/s
THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY
TANGENTIAL_COORDS = 1
# read config ini
config = cp.ConfigParser()
config_path = os.path.join('config.ini')
config.read(config_path)
sim_path = config.get('directories', 'simulation')
exp_path = config.get('directories', 'experiment')
out_path = config.get('directories', 'output')
# create output directory
out_dir = os.path.join(out_path)
if not os.path.exists(out_dir):
os.mkdir(out_dir)
h5path = os.path.join(sim_path)
h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
# read time
relative_time = np.array(h5file['relativeTime'])
real_time = relative_time * FINAL_TIME
velocityProxy<- h5file['/frictionalBoundary/velocity']
## We are interested in an enlarged time range around actual events here,
## (and no other quantities!) hence we pass a very low velocity here.
quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
indices = 1:dim(velocityProxy)[1], 1)
quakes$beginning <- realTime[quakes$beginningIndex]
quakes$ending <- realTime[quakes$endingIndex]
quakes$duration <- quakes$ending - quakes$beginning
numQuakes <- nrow(quakes)
relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
quakes[[numQuakes, 'ending']]), f=0.02)
threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
plotMask <- threeQuakeTimeMask
plotIndices<- which(plotMask)
plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
write(relaxedTime[[1]],
file.path(directories[['output']],
paste.(pasteColon('timeframe', 'min', 'threequakes',
basedir), 'tex')))
write(relaxedTime[[2]],
file.path(directories[['output']],
paste.(pasteColon('timeframe', 'max', 'threequakes',
basedir), 'tex')))
timeWindow = realTime[plotIndices]
relativeTau <- h5file['relativeTimeIncrement'][]
fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
multiGridIterations <- h5file['/iterations/multiGrid/final'][]
write.csv(data.frame(time = timeWindow,
timeIncrement = finalTime * relativeTau[plotIndices],
fixedPointIterations = fixedPointIterations[plotIndices],
multiGridIterations = multiGridIterations[plotIndices]),
file.path(directories[['output']],
paste.(pasteColon('2d-performance', basedir), 'csv')),
row.names = FALSE, quote = FALSE)
h5::h5close(h5file)
add_subdirectory("data-structures") add_subdirectory("data-structures")
add_subdirectory("factories") add_subdirectory("factories")
add_subdirectory("io") add_subdirectory("io")
add_subdirectory("multi-threading")
add_subdirectory("problem-data") add_subdirectory("problem-data")
add_subdirectory("spatial-solving") add_subdirectory("spatial-solving")
add_subdirectory("tests") add_subdirectory("tests")
......
...@@ -4,7 +4,8 @@ ...@@ -4,7 +4,8 @@
#include <dune/istl/scaledidmatrix.hh> #include <dune/istl/scaledidmatrix.hh>
#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh> #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/boundaryoperatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh> #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh> #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh> #include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
...@@ -45,11 +46,9 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass( ...@@ -45,11 +46,9 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
BoundaryPatch<GridView> const &frictionalBoundary, BoundaryPatch<GridView> const &frictionalBoundary,
ScalarMatrix &frictionalBoundaryMass) const { ScalarMatrix &frictionalBoundaryMass) const {
BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis, MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis> localMassAssembler;
LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const const BoundaryOperatorAssembler<VertexBasis, VertexBasis> boundaryAssembler(vertexBasis, vertexBasis, frictionalBoundary);
frictionalBoundaryMassAssembler(frictionalBoundary); boundaryAssembler.assemble(localMassAssembler, frictionalBoundaryMass);
vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
frictionalBoundaryMass);
} }
template <class GridView, int dimension> template <class GridView, int dimension>
......
...@@ -22,6 +22,7 @@ ...@@ -22,6 +22,7 @@
template <class GridView, int dimension> template <class GridView, int dimension>
class MyAssembler { class MyAssembler {
public: public:
using GV = GridView;
using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>; using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
using Matrix = using Matrix =
Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>; Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
......
...@@ -3,14 +3,17 @@ add_subdirectory("hdf5") ...@@ -3,14 +3,17 @@ add_subdirectory("hdf5")
add_custom_target(tectonic_dune_io SOURCES add_custom_target(tectonic_dune_io SOURCES
hdf5-bodywriter.hh hdf5-bodywriter.hh
hdf5-writer.hh hdf5-writer.hh
io-handler.hh
uniform-grid-writer.cc uniform-grid-writer.cc
vtk-bodywriter.hh
vtk.hh vtk.hh
vtk.cc
) )
#install headers #install headers
install(FILES install(FILES
hdf5-bodywriter.hh hdf5-bodywriter.hh
hdf5-levelwriter.hh hdf5-levelwriter.hh
io-handler.hh
vtk-bodywriter.hh
vtk.hh vtk.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
...@@ -18,7 +18,7 @@ class HDF5BodyWriter { ...@@ -18,7 +18,7 @@ class HDF5BodyWriter {
public: public:
using VertexCoordinates = typename ProgramState::Vector; using VertexCoordinates = typename ProgramState::Vector;
using Patch = typename FrictionalBoundaryWriter<GridView>::Patch; using Patch = typename FrictionalBoundaryWriter<GridView>::Patch;
using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >; //using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
using LocalVector = typename VertexCoordinates::block_type; using LocalVector = typename VertexCoordinates::block_type;
...@@ -29,8 +29,8 @@ class HDF5BodyWriter { ...@@ -29,8 +29,8 @@ class HDF5BodyWriter {
const size_t bodyID, const size_t bodyID,
const VertexCoordinates& vertexCoordinates, const VertexCoordinates& vertexCoordinates,
const VertexBasis& vertexBasis, const VertexBasis& vertexBasis,
const Patch& frictionPatch, const Patch& frictionPatch) :
const WeakPatches& weakPatches) : // const WeakPatches& weakPatches) :
id_(bodyID), id_(bodyID),
/*#if MY_DIM == 3 /*#if MY_DIM == 3
patchInfoWriters_(patchCount_), patchInfoWriters_(patchCount_),
......
#ifndef SRC_HDF5_WRITER_HH #ifndef SRC_HDF5_WRITER_HH
#define SRC_HDF5_WRITER_HH #define SRC_HDF5_WRITER_HH
#include <dune/fufem/hdf5/file.hh> #include <vector>
#include "hdf5/iteration-writer.hh" #include "hdf5/iteration-writer.hh"
#include "hdf5/time-writer.hh" #include "hdf5/time-writer.hh"
#include "hdf5-bodywriter.hh" #include "hdf5-bodywriter.hh"
#include <dune/fufem/hdf5/file.hh>
template <class ProgramState, class VertexBasis, class GridView> template <class ProgramState, class VertexBasis, class GridView>
class HDF5Writer { class HDF5Writer {
public: public:
...@@ -14,15 +16,15 @@ class HDF5Writer { ...@@ -14,15 +16,15 @@ class HDF5Writer {
using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>; using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
using VertexBases = std::vector<const VertexBasis*>; using VertexBases = std::vector<const VertexBasis*>;
using FrictionPatches = std::vector<const typename HDF5BodyWriter::Patch*>; using FrictionPatches = std::vector<const typename HDF5BodyWriter::Patch*>;
using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>; //using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
//friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>; //friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
HDF5Writer(HDF5::File& file, HDF5Writer(HDF5::File& file,
const VertexCoordinates& vertexCoordinates, const VertexCoordinates& vertexCoordinates,
const VertexBases& vertexBases, const VertexBases& vertexBases,
const FrictionPatches& frictionPatches, const FrictionPatches& frictionPatches)
const WeakPatches& weakPatches) //const WeakPatches& weakPatches)
: file_(file), : file_(file),
iterationWriter_(file_), iterationWriter_(file_),
timeWriter_(file_), timeWriter_(file_),
...@@ -30,7 +32,7 @@ class HDF5Writer { ...@@ -30,7 +32,7 @@ class HDF5Writer {
for (size_t i=0; i<frictionPatches_.size(); i++) { for (size_t i=0; i<frictionPatches_.size(); i++) {
if (frictionPatches_[i]->numVertices() > 0) if (frictionPatches_[i]->numVertices() > 0)
bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i], weakPatches[i])); bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i])); //, weakPatches[i]));
} }
} }
......
add_custom_target(tectonic_dune_io_hdf5 SOURCES add_custom_target(tectonic_dune_io_hdf5 SOURCES
frictionalboundary-writer.hh frictionalboundary-writer.hh
frictionalboundary-writer.cc #frictionalboundary-writer.cc
iteration-writer.hh iteration-writer.hh
iteration-writer.cc #iteration-writer.cc
patchinfo-writer.hh #patchinfo-writer.hh
patchinfo-writer.cc #patchinfo-writer.cc
restartbody-io.hh
restart-io.hh restart-io.hh
restart-io.cc
restrict.hh restrict.hh
surface-writer.hh surface-writer.hh
surface-writer.cc
time-writer.hh time-writer.hh
time-writer.cc
) )
#install headers #install headers
install(FILES install(FILES
frictionalboundary-writer.hh frictionalboundary-writer.hh
iteration-writer.hh iteration-writer.hh
patchinfo-writer.hh #patchinfo-writer.hh
restartbody-io.hh
restart-io.hh restart-io.hh
restrict.hh restrict.hh
surface-writer.hh surface-writer.hh
......
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/hdf5/singletonwriter.hh>
#include "frictionalboundary-writer.hh"
#include "restrict.hh"
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(),
Vector::block_type::dimension),
frictionalBoundaryVelocityWriter_(group_, "velocity",
frictionalBoundary.numVertices(),
Vector::block_type::dimension),
frictionalBoundaryStateWriter_(group_, "state",
frictionalBoundary.numVertices()),
frictionalBoundaryCoefficientWriter_(group_, "coefficient",
frictionalBoundary.numVertices()) {
auto const frictionalBoundaryCoordinates =
restrictToSurface(vertexCoordinates, frictionalBoundary);
HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
group_, "coordinates", frictionalBoundaryCoordinates.size(),
Vector::block_type::dimension);
setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
}
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(u, frictionalBoundary_);
addEntry(frictionalBoundaryDisplacementWriter_, timeStep, frictionalBoundaryDisplacements);
auto const frictionalBoundaryVelocities = restrictToSurface(v, frictionalBoundary_);
addEntry(frictionalBoundaryVelocityWriter_, timeStep, frictionalBoundaryVelocities);
auto const frictionalBoundaryStates = restrictToSurface(alpha, frictionalBoundary_);
addEntry(frictionalBoundaryStateWriter_, timeStep, frictionalBoundaryStates);
auto const frictionalBoundaryCoefficient = restrictToSurface(frictionCoeff, frictionalBoundary_);
addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
}
#include "frictionalboundary-writer_tmpl.cc"
#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH #ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
#define 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/file.hh>
#include <dune/fufem/hdf5/sequenceio.hh> #include <dune/fufem/hdf5/sequenceio.hh>
#include <dune/fufem/hdf5/singletonwriter.hh>
#include "restrict.hh"
template <class GridView> class FrictionalBoundaryWriter { template <class GridView> class FrictionalBoundaryWriter {
public: public:
...@@ -11,10 +14,42 @@ template <class GridView> class FrictionalBoundaryWriter { ...@@ -11,10 +14,42 @@ template <class GridView> class FrictionalBoundaryWriter {
template <class Vector> template <class Vector>
FrictionalBoundaryWriter(HDF5::Group& group, const Vector& vertexCoordinates, FrictionalBoundaryWriter(HDF5::Group& group, const Vector& vertexCoordinates,
const Patch& frictionalBoundary); const Patch& frictionalBoundary) : group_(group),
frictionalBoundary_(frictionalBoundary),
frictionalBoundaryDisplacementWriter_(group_, "displacement",
frictionalBoundary.numVertices(),
Vector::block_type::dimension),
frictionalBoundaryVelocityWriter_(group_, "velocity",
frictionalBoundary.numVertices(),
Vector::block_type::dimension),
frictionalBoundaryStateWriter_(group_, "state",
frictionalBoundary.numVertices()),
frictionalBoundaryCoefficientWriter_(group_, "coefficient",
frictionalBoundary.numVertices()) {
auto const frictionalBoundaryCoordinates =
restrictToSurface(vertexCoordinates, frictionalBoundary);
HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
group_, "coordinates", frictionalBoundaryCoordinates.size(),
Vector::block_type::dimension);
setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
}
template <class Vector, class ScalarVector> template <class Vector, class ScalarVector>
void write(const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff); void write(const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff) {
auto const frictionalBoundaryDisplacements = restrictToSurface(u, frictionalBoundary_);
addEntry(frictionalBoundaryDisplacementWriter_, timeStep, frictionalBoundaryDisplacements);
auto const frictionalBoundaryVelocities = restrictToSurface(v, frictionalBoundary_);
addEntry(frictionalBoundaryVelocityWriter_, timeStep, frictionalBoundaryVelocities);
auto const frictionalBoundaryStates = restrictToSurface(alpha, frictionalBoundary_);
addEntry(frictionalBoundaryStateWriter_, timeStep, frictionalBoundaryStates);
auto const frictionalBoundaryCoefficient = restrictToSurface(frictionCoeff, frictionalBoundary_);
addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
}
private: private:
HDF5::Group& group_; HDF5::Group& group_;
...@@ -26,4 +61,5 @@ template <class GridView> class FrictionalBoundaryWriter { ...@@ -26,4 +61,5 @@ template <class GridView> class FrictionalBoundaryWriter {
HDF5::SequenceIO<1> frictionalBoundaryStateWriter_; HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_; HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
}; };
#endif #endif
#include "../../explicitvectors.hh"
#include "../../explicitgrid.hh"
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
//#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);
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "iteration-writer.hh"
IterationWriter::IterationWriter(HDF5::Grouplike &file)
: group_(file, "iterations"),
fpiSubGroup_(group_, "fixedPoint"),
mgSubGroup_(group_, "multiGrid"),
finalMGIterationWriter_(mgSubGroup_, "final"),
finalFPIIterationWriter_(fpiSubGroup_, "final"),
totalMGIterationWriter_(mgSubGroup_, "total"),
totalFPIIterationWriter_(fpiSubGroup_, "total") {}
void IterationWriter::write(size_t timeStep,
IterationRegister const &iterationCount) {
addEntry(finalMGIterationWriter_, timeStep,
iterationCount.finalCount.multigridIterations);
addEntry(finalFPIIterationWriter_, timeStep,
iterationCount.finalCount.iterations);
addEntry(totalMGIterationWriter_, timeStep,
iterationCount.totalCount.multigridIterations);
addEntry(totalFPIIterationWriter_, timeStep,
iterationCount.totalCount.iterations);
}
#ifndef SRC_HDF_ITERATION_WRITER_HH #ifndef SRC_HDF5_ITERATION_WRITER_HH
#define SRC_HDF_ITERATION_WRITER_HH #define SRC_HDF5_ITERATION_WRITER_HH
#include <dune/fufem/hdf5/sequenceio.hh>
#include "../../time-stepping/adaptivetimestepper.hh" #include "../../time-stepping/adaptivetimestepper.hh"
#include <dune/fufem/hdf5/file.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
class IterationWriter { class IterationWriter {
public: public:
IterationWriter(HDF5::Grouplike &file); IterationWriter(HDF5::Grouplike& file)
: group_(file, "iterations"),
fpiSubGroup_(group_, "fixedPoint"),
mgSubGroup_(group_, "multiGrid"),
finalMGIterationWriter_(mgSubGroup_, "final"),
finalFPIIterationWriter_(fpiSubGroup_, "final"),
totalMGIterationWriter_(mgSubGroup_, "total"),
totalFPIIterationWriter_(fpiSubGroup_, "total") {}
void write(size_t timeStep, IterationRegister const &iterationCount); void write(size_t timeStep, const IterationRegister& iterationCount) {
addEntry(finalMGIterationWriter_, timeStep,
iterationCount.finalCount.multigridIterations);
addEntry(finalFPIIterationWriter_, timeStep,
iterationCount.finalCount.iterations);
addEntry(totalMGIterationWriter_, timeStep,
iterationCount.totalCount.multigridIterations);
addEntry(totalFPIIterationWriter_, timeStep,
iterationCount.totalCount.iterations);
}
private: private:
HDF5::Group group_; HDF5::Group group_;
...@@ -21,4 +38,5 @@ class IterationWriter { ...@@ -21,4 +38,5 @@ class IterationWriter {
HDF5::SequenceIO<0, size_t> totalMGIterationWriter_; HDF5::SequenceIO<0, size_t> totalMGIterationWriter_;
HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_; HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_;
}; };
#endif #endif
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/grid/hierarchic-approximation.hh>
#include <dune/fufem/hdf5/singletonwriter.hh>
#include "patchinfo-writer.hh"
template <class LocalVector, class GridView>
GridEvaluator<LocalVector, GridView>::GridEvaluator(const CuboidGeometry& cuboidGeometry,
ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
double const bufferWidth = 0.05 * cuboidGeometry.lengthScale();
auto const xminmax = std::minmax_element(
weakPatch.vertices.begin(), weakPatch.vertices.end(),
[](LocalVector const &a, LocalVector const &b) { return a[0] < b[0]; });
double const xmin = (*xminmax.first)[0] - bufferWidth;
double const xspan = (*xminmax.second)[0] + bufferWidth - xmin;
double const zmin = -cuboidGeometry.depth() / 2.0;
double const zspan = cuboidGeometry.depth();
double spatialResolution = 0.01 * cuboidGeometry.lengthScale();
size_t const xsteps = std::round(xspan / spatialResolution);
size_t const zsteps = std::round(zspan / spatialResolution);
xCoordinates.resize(xsteps + 1);
for (size_t xi = 0; xi <= xsteps; xi++)
xCoordinates[xi] = xmin + xi * xspan / xsteps;
zCoordinates.resize(zsteps + 1);
for (size_t zi = 0; zi <= zsteps; zi++)
zCoordinates[zi] = zmin + zi * zspan / zsteps;
HierarchicApproximation<typename GridView::Grid, GridView> const
hApproximation(gridView.grid(), gridView, 1e-6 * cuboidGeometry.lengthScale());
LocalVector global(0);
localInfo.resize(xsteps + 1);
for (size_t xi = 0; xi < xCoordinates.size(); ++xi) {
localInfo[xi].resize(zsteps + 1);
for (size_t zi = 0; zi < zCoordinates.size(); ++zi) {
global[0] = xCoordinates[xi];
global[2] = zCoordinates[zi];
auto const element = hApproximation.findEntity(global);
localInfo[xi][zi] =
std::make_pair(element, element.geometry().local(global));
}
}
}
template <class LocalVector, class GridView>
template <class Function>
Dune::Matrix<typename Function::RangeType>
GridEvaluator<LocalVector, GridView>::evaluate(Function const &f) const {
Dune::Matrix<typename Function::RangeType> ret(xCoordinates.size(),
zCoordinates.size());
for (size_t xi = 0; xi < localInfo.size(); ++xi) {
auto const &localInfoX = localInfo[xi];
for (size_t zi = 0; zi < localInfoX.size(); ++zi) {
auto const &localInfoXZ = localInfoX[zi];
f.evaluateLocal(localInfoXZ.first, localInfoXZ.second, ret[xi][zi]);
}
}
return ret;
}
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,
size_t const id)
: id_(id),
group_(file, "weakPatchGrid"+std::to_string(id_)),
vertexBasis_(vertexBasis),
gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
weakPatchGridVelocityWriter_(
group_, "velocity", gridEvaluator_.xCoordinates.size(),
gridEvaluator_.zCoordinates.size(), Vector::block_type::dimension) {
HDF5::SingletonWriter<1> weakPatchGridXCoordinateWriter(
group_, "xCoordinates", gridEvaluator_.xCoordinates.size());
setEntry(weakPatchGridXCoordinateWriter, gridEvaluator_.xCoordinates);
HDF5::SingletonWriter<1> weakPatchGridZCoordinateWriter(
group_, "zCoordinates", gridEvaluator_.zCoordinates.size());
setEntry(weakPatchGridZCoordinateWriter, gridEvaluator_.zCoordinates);
}
template <class ProgramState, class VertexBasis, class GridView>
void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
ProgramState const &programState) {
// TODO
BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]);
auto const gridVelocity = gridEvaluator_.evaluate(velocity);
addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
}
#include "patchinfo-writer_tmpl.cc"
#ifndef SRC_HDF_PATCHINFO_WRITER_HH
#define SRC_HDF_PATCHINFO_WRITER_HH
#include <dune/istl/matrix.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
#include "../../problem-data/grid/cuboidgeometry.hh"
template <class LocalVector, class GridView> class GridEvaluator {
using Element = typename GridView::Grid::template Codim<0>::Entity;
using CuboidGeometry = CuboidGeometry<typename LocalVector::field_type>;
public:
GridEvaluator(const CuboidGeometry& cuboidGeometry,
const ConvexPolyhedron<LocalVector>& weakPatch,
const GridView& gridView);
template <class Function>
Dune::Matrix<typename Function::RangeType> evaluate(Function const &f) const;
Dune::BlockVector<Dune::FieldVector<double, 1>> xCoordinates;
Dune::BlockVector<Dune::FieldVector<double, 1>> zCoordinates;
private:
std::vector<std::vector<std::pair<Element, LocalVector>>> localInfo;
};
template <class ProgramState, class VertexBasis, class GridView>
class PatchInfoWriter {
using Vector = typename ProgramState::Vector;
using LocalVector = typename Vector::block_type;
using Patch = BoundaryPatch<GridView>;
public:
PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
Patch const &frictionalBoundary,
ConvexPolyhedron<LocalVector> const &weakPatch,
size_t const id);
void write(ProgramState const &programState);
private:
const size_t id_;
HDF5::Group group_;
VertexBasis const &vertexBasis_;
GridEvaluator<LocalVector, GridView> const gridEvaluator_;
HDF5::SequenceIO<3> weakPatchGridVelocityWriter_;
};
#endif
#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>;
#ifndef SRC_HDF_RESTART_HDF_HH #ifndef SRC_IO_HDF_RESTART_IO_HH
#define SRC_HDF_RESTART_HDF_HH #define SRC_IO_HDF_RESTART_IO_HH
#include <vector>
#include <dune/fufem/hdf5/file.hh> #include <dune/fufem/hdf5/file.hh>
#include <dune/fufem/hdf5/sequenceio.hh> #include <dune/fufem/hdf5/sequenceio.hh>
template <class ProgramState> class RestartBodyIO { #include "restartbody-io.hh"
public:
RestartBodyIO(HDF5::Grouplike& group, const size_t vertexCount, const size_t id);
void write(const 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 { 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: private:
std::vector<RestartBodyIO<ProgramState>* > bodyIO_; std::vector<RestartBodyIO<ProgramState>* > bodyIO_;
HDF5::SequenceIO<0> relativeTimeWriter_; HDF5::SequenceIO<0> relativeTimeWriter_;
HDF5::SequenceIO<0> relativeTimeIncrementWriter_; HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
public:
RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts)
: bodyIO_(vertexCounts.size()),
relativeTimeWriter_(group, "relativeTime"),
relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {
for (size_t i=0; i<bodyIO_.size(); i++) {
bodyIO_[i] = new RestartBodyIO<ProgramState>(group, vertexCounts[i], i);
}
}
~RestartIO() {
for (size_t i=0; i<bodyIO_.size(); i++) {
delete bodyIO_[i];
}
}
void write(ProgramState const &programState) {
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,
programState.relativeTau);
}
void read(size_t timeStep, ProgramState& programState) {
assert(programState.size() == bodyIO_.size());
programState.timeStep = timeStep;
for (size_t i=0; i<bodyIO_.size(); i++) {
bodyIO_[i]->read(timeStep, programState);
}
readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
}
}; };
#endif #endif
#include "../../explicitvectors.hh"
#include "../../data-structures/program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>;
template class RestartBodyIO<MyProgramState>;
template class RestartIO<MyProgramState>;
#ifdef HAVE_CONFIG_H #ifndef SRC_HDF_RESTART_BODY_HDF_HH
#include "config.h" #define SRC_HDF_RESTART_BODY_HDF_HH
#endif
#include <vector>
#include <dune/fufem/hdf5/file.hh>
#include <dune/fufem/hdf5/sequenceio.hh>
#include "restart-io.hh" template <class ProgramState> class RestartBodyIO {
private:
const size_t id_;
template <class ProgramState> HDF5::SequenceIO<2> displacementWriter_;
RestartBodyIO<ProgramState>::RestartBodyIO(HDF5::Grouplike &group, const size_t vertexCount, const size_t id) HDF5::SequenceIO<2> velocityWriter_;
HDF5::SequenceIO<2> accelerationWriter_;
HDF5::SequenceIO<1> stateWriter_;
HDF5::SequenceIO<1> weightedNormalStressWriter_;
public:
RestartBodyIO(HDF5::Grouplike &group, const size_t vertexCount, const size_t id)
: id_(id), : id_(id),
displacementWriter_(group, "displacement"+std::to_string(id_), vertexCount, displacementWriter_(group, "displacement"+std::to_string(id_), vertexCount,
ProgramState::Vector::block_type::dimension), ProgramState::Vector::block_type::dimension),
...@@ -16,9 +28,8 @@ RestartBodyIO<ProgramState>::RestartBodyIO(HDF5::Grouplike &group, const size_t ...@@ -16,9 +28,8 @@ RestartBodyIO<ProgramState>::RestartBodyIO(HDF5::Grouplike &group, const size_t
stateWriter_(group, "state"+std::to_string(id_), vertexCount), stateWriter_(group, "state"+std::to_string(id_), vertexCount),
weightedNormalStressWriter_(group, "weightedNormalStress"+std::to_string(id_), vertexCount) {} weightedNormalStressWriter_(group, "weightedNormalStress"+std::to_string(id_), vertexCount) {}
template <class ProgramState>
void RestartBodyIO<ProgramState>::write(const ProgramState& programState) {
void write(const ProgramState& programState) {
addEntry(displacementWriter_, programState.timeStep, programState.u[id_]); addEntry(displacementWriter_, programState.timeStep, programState.u[id_]);
addEntry(velocityWriter_, programState.timeStep, programState.v[id_]); addEntry(velocityWriter_, programState.timeStep, programState.v[id_]);
addEntry(accelerationWriter_, programState.timeStep, programState.a[id_]); addEntry(accelerationWriter_, programState.timeStep, programState.a[id_]);
...@@ -27,10 +38,8 @@ void RestartBodyIO<ProgramState>::write(const ProgramState& programState) { ...@@ -27,10 +38,8 @@ void RestartBodyIO<ProgramState>::write(const ProgramState& programState) {
programState.weightedNormalStress[id_]); programState.weightedNormalStress[id_]);
} }
template <class ProgramState>
void RestartBodyIO<ProgramState>::read(size_t timeStep,
ProgramState& programState) {
void read(size_t timeStep, ProgramState& programState) {
readEntry(displacementWriter_, timeStep, programState.u[id_]); readEntry(displacementWriter_, timeStep, programState.u[id_]);
readEntry(velocityWriter_, timeStep, programState.v[id_]); readEntry(velocityWriter_, timeStep, programState.v[id_]);
readEntry(accelerationWriter_, timeStep, programState.a[id_]); readEntry(accelerationWriter_, timeStep, programState.a[id_]);
...@@ -39,51 +48,5 @@ void RestartBodyIO<ProgramState>::read(size_t timeStep, ...@@ -39,51 +48,5 @@ void RestartBodyIO<ProgramState>::read(size_t timeStep,
programState.weightedNormalStress[id_]); programState.weightedNormalStress[id_]);
} }
template <class ProgramState> };
RestartIO<ProgramState>::RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts) #endif
: bodyIO_(vertexCounts.size()),
relativeTimeWriter_(group, "relativeTime"),
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) {
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,
programState.relativeTau);
}
template <class ProgramState>
void RestartIO<ProgramState>::read(size_t timeStep,
ProgramState& programState) {
assert(programState.size() == bodyIO_.size());
programState.timeStep = timeStep;
for (size_t i=0; i<bodyIO_.size(); i++) {
bodyIO_[i]->read(timeStep, programState);
}
readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
}
#include "restart-io_tmpl.cc"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "restrict.hh"
#include "surface-writer.hh"
template <class ProgramState, class GridView>
SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface, size_t const id)
: id_(id),
group_(file, "surface"+std::to_string(id_)),
surface_(surface),
surfaceDisplacementWriter_(group_, "displacement", surface.numVertices(),
Vector::block_type::dimension),
surfaceVelocityWriter_(group_, "velocity", surface.numVertices(),
Vector::block_type::dimension) {
auto const surfaceCoordinates = restrictToSurface(vertexCoordinates, surface);
HDF5::SingletonWriter<2> surfaceCoordinateWriter(group_, "coordinates",
surfaceCoordinates.size(),
Vector::block_type::dimension);
setEntry(surfaceCoordinateWriter, surfaceCoordinates);
}
template <class ProgramState, class GridView>
void SurfaceWriter<ProgramState, GridView>::write(
ProgramState const &programState) {
auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_);
addEntry(surfaceDisplacementWriter_, programState.timeStep,
surfaceDisplacements);
auto const surfaceVelocities = restrictToSurface(programState.v[id_], surface_);
addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
}
#include "surface-writer_tmpl.cc"
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment