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 459 additions and 164 deletions
// Copyright Carsten Graeser 2012
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
template <class Enum> struct StringToEnum : public Dune::NotImplemented {};
template <class Enum>
typename Dune::enable_if<
!Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value,
std::istream &>::type
operator>>(std::istream &lhs, Enum &e) {
std::string s;
lhs >> s;
try {
e = StringToEnum<Enum>::convert(s);
}
catch (typename Dune::Exception) {
lhs.setstate(std::ios_base::failbit);
}
return lhs;
}
#include <dune/common/exceptions.hh>
#include "enums.hh"
template <> struct StringToEnum<Config::scheme> {
static Config::scheme convert(std::string const &s) {
if (s == "newmark")
return Config::Newmark;
if (s == "backwardEuler")
return Config::BackwardEuler;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
};
#include <dune/common/exceptions.hh>
#include "enums.hh"
template <> struct StringToEnum<Config::stateModel> {
static Config::stateModel convert(std::string const &s) {
if (s == "AgeingLaw")
return Config::AgeingLaw;
if (s == "SlipLaw")
return Config::SlipLaw;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
};
#include <dune/solvers/common/numproc.hh> // Solver::VerbosityMode
#include <dune/common/exceptions.hh>
template <> struct StringToEnum<Solver::VerbosityMode> {
static Solver::VerbosityMode convert(std::string const &s) {
if (s == "full")
return Solver::FULL;
if (s == "reduced")
return Solver::REDUCED;
if (s == "quiet")
return Solver::QUIET;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
};
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
// Copyright Carsten Graeser 2012
#include <dune/common/exceptions.hh>
#include "enumparser.hh"
template <class Enum>
typename std::enable_if<
!std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
std::istream &>::type
operator>>(std::istream &lhs, Enum &e) {
std::string s;
lhs >> s;
try {
e = StringToEnum<Enum>::convert(s);
} catch (typename Dune::Exception) {
lhs.setstate(std::ios_base::failbit);
}
return lhs;
}
Config::FrictionModel StringToEnum<Config::FrictionModel>::convert(
std::string const &s) {
if (s == "Truncated")
return Config::Truncated;
if (s == "Regularised")
return Config::Regularised;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
Config::scheme StringToEnum<Config::scheme>::convert(std::string const &s) {
if (s == "newmark")
return Config::Newmark;
if (s == "backwardEuler")
return Config::BackwardEuler;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
Config::stateModel StringToEnum<Config::stateModel>::convert(
std::string const &s) {
if (s == "AgeingLaw")
return Config::AgeingLaw;
if (s == "SlipLaw")
return Config::SlipLaw;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
Config::PatchType StringToEnum<Config::PatchType>::convert(
std::string const &s) {
if (s == "Rectangular")
return Config::Rectangular;
if (s == "Trapezoidal")
return Config::Trapezoidal;
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
template std::istream &operator>>(std::istream &lhs, Config::FrictionModel &);
template std::istream &operator>>(std::istream &lhs, Config::stateModel &);
template std::istream &operator>>(std::istream &lhs, Config::scheme &);
template std::istream &operator>>(std::istream &lhs, Config::PatchType &);
#ifndef SRC_ENUMPARSER_HH
#define SRC_ENUMPARSER_HH
// Copyright Carsten Graeser 2012
#include <type_traits>
#include <dune/solvers/solvers/solver.hh>
#include "enums.hh"
template <class Enum> struct StringToEnum : public Dune::NotImplemented {};
template <> struct StringToEnum<Config::FrictionModel> {
static Config::FrictionModel convert(std::string const &s);
};
template <> struct StringToEnum<Config::stateModel> {
static Config::stateModel convert(std::string const &s);
};
template <> struct StringToEnum<Config::scheme> {
static Config::scheme convert(std::string const &s);
};
template <> struct StringToEnum<Config::PatchType> {
static Config::PatchType convert(std::string const &s);
};
template <class Enum>
typename std::enable_if<
!std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
std::istream &>::type
operator>>(std::istream &lhs, Enum &e);
#endif
......@@ -2,18 +2,10 @@
#define SRC_ENUMS_HH
struct Config {
enum FrictionModel {
Truncated,
Regularised
};
enum stateModel {
AgeingLaw,
SlipLaw
};
enum scheme {
Newmark,
BackwardEuler
};
enum FrictionModel { Truncated, Regularised };
enum stateModel { AgeingLaw, SlipLaw };
enum scheme { Newmark, BackwardEuler };
enum PatchType { Rectangular, Trapezoidal };
};
#endif
#ifndef SRC_EXPLICITGRID_HH
#define SRC_EXPLICITGRID_HH
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wignored-qualifiers"
#include <dune/grid/alugrid.hh>
#pragma clang diagnostic pop
#include "gridselector.hh"
using Grid = Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming>;
using GridView = Grid::LeafGridView;
#endif
......@@ -6,9 +6,10 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
using LocalVector = Dune::FieldVector<double, DIM>;
using LocalMatrix = Dune::FieldMatrix<double, DIM, DIM>;
using LocalVector = Dune::FieldVector<double, MY_DIM>;
using LocalMatrix = Dune::FieldMatrix<double, MY_DIM, MY_DIM>;
using Vector = Dune::BlockVector<LocalVector>;
using Matrix = Dune::BCRSMatrix<LocalMatrix>;
using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
#endif
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "friction_writer.hh"
#include "tobool.hh"
template <class ScalarVector, class Vector>
FrictionWriter<ScalarVector, Vector>::FrictionWriter(
Vector const &vertexCoordinates,
Dune::BitSetVector<1> const &_boundaryNodes, std::string prefix,
typename BW::Projector projector)
: BW(vertexCoordinates, _boundaryNodes, prefix, projector),
coefficientWriter(prefix + "Coefficients", std::fstream::out),
stateWriter(prefix + "Alpha", std::fstream::out) {}
template <class ScalarVector, class Vector>
FrictionWriter<ScalarVector, Vector>::~FrictionWriter() {
stateWriter.close();
coefficientWriter.close();
}
template <class ScalarVector, class Vector>
void FrictionWriter<ScalarVector, Vector>::writeOther(
ScalarVector const &coefficient, ScalarVector const &alpha) {
for (size_t i = 0; i < boundaryNodes.size(); ++i) {
if (!toBool(boundaryNodes[i]))
continue;
coefficientWriter << coefficient[i] << " ";
stateWriter << alpha[i] << " ";
}
stateWriter << std::endl;
coefficientWriter << std::endl;
}
#include "friction_writer_tmpl.cc"
#ifndef SRC_FRICTION_WRITER_HH
#define SRC_FRICTION_WRITER_HH
#include <fstream>
#include <string>
#include <dune/common/bitsetvector.hh>
#include "boundary_writer.hh"
template <class ScalarVector, class Vector>
class FrictionWriter : public BoundaryWriter<ScalarVector, Vector> {
using BW = BoundaryWriter<ScalarVector, Vector>;
public:
FrictionWriter(Vector const &vertexCoordinates,
Dune::BitSetVector<1> const &_boundaryNodes,
std::string prefix, typename BW::Projector projector);
~FrictionWriter();
void writeOther(ScalarVector const &coefficient, ScalarVector const &alpha);
private:
std::fstream coefficientWriter;
std::fstream stateWriter;
using BW::boundaryNodes;
};
#endif
#ifndef DIM
#error DIM unset
#endif
#include "explicitvectors.hh"
template class FrictionWriter<ScalarVector, Vector>;
#ifndef MY_DIM
#error MY_DIM unset
#endif
#define WANT_ALUGRID 0
#define WANT_UG 1
#define WANT_GRID WANT_UG
#if WANT_GRID == WANT_ALUGRID
#if !HAVE_ALUGRID
#error ALUGRID was requested but not found
#endif
#include <dune/grid/alugrid.hh>
using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
#elif WANT_GRID == WANT_UG
#if !HAVE_UG
#error UG was requested but not found
#endif
#include <dune/grid/uggrid.hh>
using Grid = Dune::UGGrid<MY_DIM>;
#else
#error requested a grid that does not exist
#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
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/hdf5/singletonwriter.hh>
#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)
: group_(file, "frictionalBoundary"),
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 ProgramState, class GridView>
template <class Friction>
void FrictionalBoundaryWriter<ProgramState, GridView>::write(
ProgramState const &programState, Friction &friction) {
auto const frictionalBoundaryDisplacements =
restrictToSurface(programState.u, frictionalBoundary_);
addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
frictionalBoundaryDisplacements);
auto const frictionalBoundaryVelocities =
restrictToSurface(programState.v, frictionalBoundary_);
addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
frictionalBoundaryVelocities);
auto const frictionalBoundaryStates =
restrictToSurface(programState.alpha, frictionalBoundary_);
addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
frictionalBoundaryStates);
friction.updateAlpha(programState.alpha);
auto const c = friction.coefficientOfFriction(programState.v);
auto const frictionalBoundaryCoefficient =
restrictToSurface(c, frictionalBoundary_);
addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
frictionalBoundaryCoefficient);
}
#include "frictionalboundary-writer_tmpl.cc"
#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_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;
using Patch = BoundaryPatch<GridView>;
public:
FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary);
template <class Friction>
void write(ProgramState const &programState, Friction &friction);
private:
HDF5::Group group_;
Patch const &frictionalBoundary_;
HDF5::SequenceIO<2> frictionalBoundaryDisplacementWriter_;
HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
};
#endif
#include "../explicitvectors.hh"
#include "../explicitgrid.hh"
#include "../program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>;
using MyFriction = GlobalFriction<Matrix, Vector>;
template class FrictionalBoundaryWriter<MyProgramState, GridView>;
template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
MyProgramState const &programState, MyFriction &friction);
#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
#define SRC_HDF_ITERATION_WRITER_HH
#include <dune/fufem/hdf5/sequenceio.hh>
#include "../time-stepping/adaptivetimestepper.hh"
class IterationWriter {
public:
IterationWriter(HDF5::Grouplike &file);
void write(size_t timeStep, IterationRegister const &iterationCount);
private:
HDF5::Group group_;
HDF5::Group fpiSubGroup_;
HDF5::Group mgSubGroup_;
HDF5::SequenceIO<0, size_t> finalMGIterationWriter_;
HDF5::SequenceIO<0, size_t> finalFPIIterationWriter_;
HDF5::SequenceIO<0, size_t> totalMGIterationWriter_;
HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_;
};
#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(
ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
double const bufferWidth = 0.05 * MyGeometry::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 = -MyGeometry::depth / 2.0;
double const zspan = MyGeometry::depth;
double spatialResolution = 0.01 * MyGeometry::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 * MyGeometry::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)
: group_(file, "weakPatchGrid"),
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) {
BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
auto const gridVelocity = gridEvaluator_.evaluate(velocity);
addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
}
#include "patchinfo-writer_tmpl.cc"