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
Hide whitespace changes
Inline Side-by-side
Showing
with 791 additions and 559 deletions
dune/tectonic/utils/maxeigenvalue.hh 0 → 100644
View file @ cdeb7bf9
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
/** \file
* \brief max eigenvalue computations for the FieldMatrix class
*/
#include <iostream>
#include <cmath>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
return matrix[0][0];
}
/** \brief calculates the eigenvalues of a symmetric field matrix
\param[in] matrix matrix eigenvalues are calculated for
\param[out] max eigenvalue
*/
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
const K p2 = p - matrix[1][1];
K q = p2 * p2 + matrix[1][0] * matrix[0][1];
if( q < 0 && q > -1e-14 ) q = 0;
if (q < 0) {
std::cout << matrix << std::endl;
// Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
}
// get square root
q = std::sqrt(q);
return p + q;
}
/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
\param[in] matrix eigenvalues are calculated for
\param[out] max eigenvalue
\note If the input matrix is not symmetric the behavior of this method is undefined.
*/
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
Dune::FieldVector<K, 3> eigenvalues;
Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
return eigenvalues[0];
}
/** @} end documentation */
#endif
dune/tectonic/utils/reductionfactors.hh 0 → 100644
View file @ cdeb7bf9
extern std::vector<std::vector<double>> allReductionFactors;
src/tobool.hh → dune/tectonic/utils/tobool.hh
View file @ cdeb7bf9
......@@ -4,7 +4,7 @@
#include <dune/common/bitsetvector.hh>
template <class Alloc>
bool toBool(Dune::BitSetVectorConstReference<1, Alloc> x) {
bool toBool(const Dune::BitSetVectorConstReference<1, Alloc> x) {
return x[0];
}
#endif
src/CMakeLists.txt
View file @ cdeb7bf9
set(SW_SOURCE_FILES
assemblers.cc
enumparser.cc
hdf5/frictionalboundary-writer.cc
hdf5/iteration-writer.cc
hdf5/patchinfo-writer.cc
hdf5/restart-io.cc
hdf5/surface-writer.cc
hdf5/time-writer.cc
one-body-problem-data/mygeometry.cc
one-body-problem-data/mygrid.cc
one-body-problem.cc
spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc
time-stepping/adaptivetimestepper.cc
time-stepping/coupledtimestepper.cc
time-stepping/rate.cc
time-stepping/rate/rateupdater.cc
time-stepping/state.cc
vtk.cc
)
add_subdirectory("foam")
add_subdirectory("multi-body-problem")
add_subdirectory("strikeslip")
set(UGW_SOURCE_FILES
assemblers.cc # FIXME
one-body-problem-data/mygrid.cc
uniform-grid-writer.cc
vtk.cc
../dune/tectonic/assemblers.cc # FIXME
#../dune/tectonic/io/uniform-grid-writer.cc
../dune/tectonic/io/vtk.cc
../dune/tectonic/problem-data/grid/mygrids.cc
)
foreach(_dim 2 3)
set(_sw_target one-body-problem-${_dim}D)
set(_ugw_target uniform-grid-writer-${_dim}D)
add_executable(${_sw_target} ${SW_SOURCE_FILES})
add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
add_dune_ug_flags(${_sw_target})
add_dune_ug_flags(${_ugw_target})
#set(_ugw_target uniform-grid-writer-${_dim}D)
add_dune_hdf5_flags(${_sw_target})
#add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
#add_dune_ug_flags(${_ugw_target})
set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
#set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
endforeach()
src/assemblers.hh deleted 100644 → 0
View file @ 8fe950bc
#ifndef SRC_ASSEMBLERS_HH
#define SRC_ASSEMBLERS_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/assemblers/assembler.hh>
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wsign-compare"
#include <dune/fufem/functionspacebases/p0basis.hh>
#pragma clang diagnostic pop
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/globalfrictiondata.hh>
#include "enums.hh"
template <class GridView, int dimension> class MyAssembler {
public:
using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
using Matrix =
Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
using CellBasis = P0Basis<GridView, double>;
using VertexBasis = P1NodalBasis<GridView, double>;
CellBasis const cellBasis;
VertexBasis const vertexBasis;
GridView const &gridView;
private:
using Grid = typename GridView::Grid;
using LocalVector = typename Vector::block_type;
using LocalScalarVector = typename ScalarVector::block_type;
using LocalCellBasis = typename CellBasis::LocalFiniteElement;
using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
Assembler<CellBasis, CellBasis> cellAssembler;
Assembler<VertexBasis, VertexBasis> vertexAssembler;
public:
MyAssembler(GridView const &gridView);
void assembleFrictionalBoundaryMass(
BoundaryPatch<GridView> const &frictionalBoundary,
ScalarMatrix &frictionalBoundaryMass) const;
void assembleMass(Dune::VirtualFunction<
LocalVector, LocalScalarVector> const &densityFunction,
Matrix &M) const;
void assembleElasticity(double E, double nu, Matrix &A) const;
void assembleViscosity(
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
shearViscosity,
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
bulkViscosity,
Matrix &C) const;
void assembleBodyForce(
Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
Vector &f) const;
void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
Vector &f,
Dune::VirtualFunction<double, double> const &neumann,
double relativeTime) const;
void assembleWeightedNormalStress(
BoundaryPatch<GridView> const &frictionalBoundary,
ScalarVector &weightedNormalStress, double youngModulus,
double poissonRatio, Vector const &displacement) const;
std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
Config::FrictionModel frictionModel,
BoundaryPatch<GridView> const &frictionalBoundary,
GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &weightedNormalStress) const;
void assembleVonMisesStress(double youngModulus, double poissonRatio,
Vector const &u, ScalarVector &stress) const;
};
#endif
src/explicitgrid.hh deleted 100644 → 0
View file @ 8fe950bc
#ifndef SRC_EXPLICITGRID_HH
#define SRC_EXPLICITGRID_HH
#include "gridselector.hh"
using GridView = Grid::LeafGridView;
#endif
src/foam/CMakeLists.txt 0 → 100644
View file @ cdeb7bf9
add_custom_target(tectonic_src_foam SOURCES
foam.cfg
foam-2D.cfg
)
set(FOAM_SOURCE_FILES
../../dune/tectonic/assemblers.cc
../../dune/tectonic/data-structures/body/body.cc
../../dune/tectonic/data-structures/network/levelcontactnetwork.cc
../../dune/tectonic/data-structures/network/contactnetwork.cc
../../dune/tectonic/data-structures/enumparser.cc
../../dune/tectonic/factories/twoblocksfactory.cc
#../../dune/tectonic/io/vtk.cc
#../../dune/tectonic/io/hdf5/frictionalboundary-writer.cc
#../../dune/tectonic/io/hdf5/iteration-writer.cc
#../../dune/tectonic/io/hdf5/patchinfo-writer.cc
#../../dune/tectonic/io/hdf5/restart-io.cc
#../../dune/tectonic/io/hdf5/surface-writer.cc
#../../dune/tectonic/io/hdf5/time-writer.cc
../../dune/tectonic/problem-data/grid/cuboidgeometry.cc
../../dune/tectonic/problem-data/grid/mygrids.cc
../../dune/tectonic/problem-data/grid/simplexmanager.cc
../../dune/tectonic/spatial-solving/solverfactory.cc
../../dune/tectonic/spatial-solving/fixedpointiterator.cc
../../dune/tectonic/time-stepping/coupledtimestepper.cc
../../dune/tectonic/time-stepping/adaptivetimestepper.cc
../../dune/tectonic/time-stepping/rate.cc
../../dune/tectonic/time-stepping/rate/rateupdater.cc
../../dune/tectonic/time-stepping/state.cc
../../dune/tectonic/time-stepping/uniformtimestepper.cc
foam.cc
)
foreach(_dim 2 3)
set(_foam_target foam-${_dim}D)
add_executable(${_foam_target} ${FOAM_SOURCE_FILES})
add_dune_ug_flags(${_foam_target})
add_dune_hdf5_flags(${_foam_target})
set_property(TARGET ${_foam_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
#set_property(TARGET ${_foam_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
endforeach()
src/foam/foam-2D.cfg 0 → 100644
View file @ cdeb7bf9
# -*- mode:conf -*-
[body0]
smallestDiameter = 6.25e-2 # 2e-3 [m]
[body1]
smallestDiameter = 6.25e-2 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
[u0.solver]
tolerance = 1e-8
[a0.solver]
tolerance = 1e-8
[v.solver]
tolerance = 1e-8
[v.fpi]
tolerance = 1e-10
[solver.tnnmg.preconditioner.basesolver]
tolerance = 1e-10
[solver.tnnmg.preconditioner.patchsolver]
tolerance = 1e-10
src/foam/foam.cc 0 → 100644
View file @ cdeb7bf9
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_IPOPT
#undef HAVE_IPOPT
#endif
#include <atomic>
#include <cmath>
#include <csignal>
#include <exception>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <memory>
#include <dune/common/bitsetvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/formatstring.hh>
#include <dune/fufem/hdf5/file.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/contact/common/deformedcontinuacomplex.hh>
#include <dune/contact/common/couplingpair.hh>
#include <dune/contact/projections/normalprojection.hh>
#include <dune/tnnmg/functionals/quadraticfunctional.hh>
#include <dune/tectonic/assemblers.hh>
#include <dune/tectonic/gridselector.hh>
#include <dune/tectonic/explicitgrid.hh>
#include <dune/tectonic/explicitvectors.hh>
#include <dune/tectonic/data-structures/enumparser.hh>
#include <dune/tectonic/data-structures/enums.hh>
#include <dune/tectonic/data-structures/network/contactnetwork.hh>
#include <dune/tectonic/data-structures/matrices.hh>
#include <dune/tectonic/data-structures/program_state.hh>
#include <dune/tectonic/data-structures/friction/globalfriction.hh>
#include <dune/tectonic/factories/twoblocksfactory.hh>
#include <dune/tectonic/io/io-handler.hh>
#include <dune/tectonic/io/hdf5-writer.hh>
#include <dune/tectonic/io/hdf5/restart-io.hh>
#include <dune/tectonic/io/vtk.hh>
#include <dune/tectonic/problem-data/bc.hh>
#include <dune/tectonic/problem-data/mybody.hh>
#include <dune/tectonic/problem-data/grid/mygrids.hh>
#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
//#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
#include <dune/tectonic/spatial-solving/solverfactory.hh>
#include <dune/tectonic/spatial-solving/makelinearsolver.hh>
#include <dune/tectonic/time-stepping/uniformtimestepper.hh>
#include <dune/tectonic/time-stepping/rate.hh>
#include <dune/tectonic/time-stepping/state.hh>
#include <dune/tectonic/time-stepping/updaters.hh>
#include <dune/tectonic/utils/debugutils.hh>
#include <dune/tectonic/utils/diameter.hh>
#include <dune/tectonic/utils/geocoordinate.hh>
#include <dune/matrix-vector/axpy.hh>
#include <dune/contact/assemblers/nbodyassembler.hh>
#include "dune/tectonic/spatial-solving/tnnmg/zerononlinearity.hh"
// for getcwd
#include <unistd.h>
//#include <tbb/tbb.h> //TODO multi threading preconditioner?
//#include <pthread.h>
#include <dune/tectonic/utils/reductionfactors.hh>
std::vector<std::vector<double>> allReductionFactors;
size_t const dims = MY_DIM;
const std::string sourcePath = "/home/joscha/software/dune/dune-tectonic/src/foam/";
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree(sourcePath + "foam.cfg", parset);
Dune::ParameterTreeParser::readINITree(
Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
return parset;
}
static std::atomic<bool> terminationRequested(false);
void handleSignal(int signum) { terminationRequested = true; }
int main(int argc, char *argv[]) {
using BlocksFactory = TwoBlocksFactory<Grid, Vector>;
using ContactNetwork = typename BlocksFactory::ContactNetwork;
using MyProgramState = ProgramState<Vector, ScalarVector>;
using Assembler = MyAssembler<DefLeafGridView, dims>;
using IOHandler = IOHandler<Assembler, ContactNetwork, MyProgramState>;
std::unique_ptr<IOHandler> ioHandler;
try {
Dune::MPIHelper::instance(argc, argv);
char buffer[256];
char *val = getcwd(buffer, sizeof(buffer));
if (val) {
std::cout << buffer << std::endl;
std::cout << argv[0] << std::endl;
}
auto const parset = getParameters(argc, argv);
auto outPath = std::filesystem::current_path();
outPath += "/output/" + parset.get<std::string>("outPath");
if (!std::filesystem::is_directory(outPath))
std::filesystem::create_directories(outPath);
const auto copyOptions = std::filesystem::copy_options::overwrite_existing;
std::filesystem::copy(sourcePath + "foam.cfg", outPath, copyOptions);
std::filesystem::copy(Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), outPath, copyOptions);
std::filesystem::current_path(outPath);
std::ofstream out("foam.log");
std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
using field_type = Matrix::field_type;
// ----------------------
// set up contact network
// ----------------------
BlocksFactory blocksFactory(parset);
blocksFactory.build();
ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
const size_t bodyCount = contactNetwork.nBodies();
for (size_t i=0; i<contactNetwork.nLevels(); i++) {
// printDofLocation(contactNetwork.body(i)->gridView());
//Vector def(contactNetwork.deformedGrids()[i]->size(dims));
//def = 1;
//deformedGridComplex.setDeformation(def, i);
const auto& level = *contactNetwork.level(i);
//for (size_t j=0; j<level.nBodies(); j++) {
//writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
//}
}
//for (size_t i=0; i<bodyCount; i++) {
//writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
//}
// ----------------------------
// assemble contactNetwork
// ----------------------------
contactNetwork.assemble();
//printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
// -----------------
// init input/output
// -----------------
std::vector<size_t> nVertices(bodyCount);
for (size_t i=0; i<bodyCount; i++) {
nVertices[i] = contactNetwork.body(i)->nVertices();
}
print(nVertices, "#dofs: ");
MyProgramState programState(nVertices);
ioHandler = std::make_unique<IOHandler>(parset.sub("io"), contactNetwork);
bool restartRead = ioHandler->read(programState);
if (!restartRead) {
programState.setupInitialConditions(parset, contactNetwork);
}
const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
/* auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
auto nonlinearity = ZeroNonlinearity();
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
Dune::ParameterTree const &_localParset) {
Vector totalX;
nBodyAssembler.nodalToTransformed(_x, totalX);
FunctionalFactory<std::decay_t<decltype(nBodyAssembler)>, decltype(nonlinearity), Matrix, Vector> functionalFactory(nBodyAssembler, nonlinearity, _matrices, _rhs);
functionalFactory.build();
auto functional = functionalFactory.functional();
NonlinearSolver<std::remove_reference_t<decltype(*functional)>, BitVector> solver(parset.sub("solver.tnnmg"), linearSolver, functional, _dirichletNodes);
solver.solve(_localParset, totalX);
nBodyAssembler.postprocess(totalX, _x);
};
programState.timeStep = parset.get<size_t>("initialTime.timeStep");
programState.relativeTime = parset.get<double>("initialTime.relativeTime");
programState.relativeTau = parset.get<double>("initialTime.relativeTau");
std::vector<Vector> ell0(bodyCount);
for (size_t i=0; i<bodyCount; i++) {
programState.u[i] = 0.0;
programState.v[i] = 0.0;
programState.a[i] = 0.0;
ell0[i].resize(programState.u[i].size());
ell0[i] = 0.0;
contactNetwork.body(i)->externalForce()(programState.relativeTime, ell0[i]);
}
*/
// Initial displacement: Start from a situation of minimal stress,
// which is automatically attained in the case [v = 0 = a].
// Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
BitVector totalDirichletNodes;
contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
size_t dof = 0;
for (size_t i=0; i<bodyCount; i++) {
const auto& body = *contactNetwork.body(i);
if (body.data()->getYoungModulus() == 0.0) {
for (; dof<body.nVertices(); dof++) {
totalDirichletNodes[dof] = true;
}
} else {
dof += body.nVertices();
}
}
using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
BoundaryNodes dirichletNodes;
contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
std::vector<const Dune::BitSetVector<1>*> frictionNodes;
contactNetwork.frictionNodes(frictionNodes);
/*
std::cout << "solving linear problem for u..." << std::endl;
{
int nVertices = contactNetwork.body(1)->nVertices();
BitVector uDirichletNodes(nVertices, false);
int idx=0;
const auto& body = contactNetwork.body(1);
using LeafBody = typename ContactNetwork::LeafBody;
std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
body->boundaryConditions("dirichlet", boundaryConditions);
if (boundaryConditions.size()>0) {
const int idxBackup = idx;
for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
const auto& nodes = boundaryConditions[bc]->boundaryNodes();
for (size_t j=0; j<nodes->size(); j++, idx++)
for (int k=0; k<dims; k++)
uDirichletNodes[idx][k] = uDirichletNodes[idx][k] or (*nodes)[j][k];
idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
}
}
for (auto i=0; i<nVertices; i++) {
if ((*frictionNodes[1])[i][0]) {
uDirichletNodes[i][1] = true;
}
}
using Norm = EnergyNorm<Matrix, Vector>;
using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
// transfer operators need to be recomputed on change due to setDeformation()
using TransferOperator = CompressedMultigridTransfer<Vector>;
using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
const auto& grid = contactNetwork.body(1)->grid();
TransferOperators transfer(grid->maxLevel());
for (size_t i = 0; i < transfer.size(); ++i)
{
// create transfer operators from level i to i+1 (note that this will only work for either uniformly refined grids or adaptive grids with RefinementType=COPY)
auto t = std::make_shared<TransferOperator>();
t->setup(*grid, i, i+1);
transfer[i] = t;
}
auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >(
*contactNetwork.matrices().elasticity[1],
programState.u[1],
ell0[1]);
linearMultigridStep->setIgnore(uDirichletNodes);
linearMultigridStep->setMGType(1, 3, 3);
linearMultigridStep->setSmoother(TruncatedBlockGSStep<Matrix, Vector>());
linearMultigridStep->setTransferOperators(transfer);
const auto& solverParset = parset.sub("u0.solver");
Dune::Solvers::LoopSolver<Vector> solver(linearMultigridStep, solverParset.get<size_t>("maximumIterations"),
solverParset.get<double>("tolerance"), Norm(*linearMultigridStep),
solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
solver.preprocess();
solver.solve();
}
//print(u, "initial u:");
// 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++) {
const auto& body = contactNetwork.body(i);
Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, programState.u[i]);
}
std::cout << "solving linear problem for a..." << std::endl;
BitVector noNodes(totalDirichletNodes.size(), false);
solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
parset.sub("a0.solver"));
// setup initial conditions in program state
programState.alpha[0] = 0;
programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha");
for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
const auto& coupling = contactNetwork.coupling(i);
const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
const auto nonmortarIdx = coupling->gridIdx_[0];
const auto& body = contactNetwork.body(nonmortarIdx);
ScalarVector frictionBoundaryStress(programState.weightedNormalStress[nonmortarIdx].size());
body->assembler()->assembleWeightedNormalStress(
contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
body->data()->getPoissonRatio(), programState.u[nonmortarIdx]);
programState.weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
} */
for (size_t i=0; i<bodyCount; i++) {
contactNetwork.body(i)->setDeformation(programState.u[i]);
}
nBodyAssembler.assembleTransferOperator();
nBodyAssembler.assembleObstacle();
// ------------------------
// assemble global friction
// ------------------------
//print(programState.weightedNormalStress, "weightedNormalStress");
contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
auto& globalFriction = contactNetwork.globalFriction();
globalFriction.updateAlpha(programState.alpha);
IterationRegister iterationCount;
ioHandler->write(programState, contactNetwork, globalFriction, iterationCount, true);
// -------------------
// Set up TNNMG solver
// -------------------
//using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
using NonlinearFactory = SolverFactory<Functional, BitVector>;
using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
StateUpdater<ScalarVector, Vector>>;
BoundaryFunctions velocityDirichletFunctions;
contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
/*for (size_t i=0; i<dirichletNodes.size(); i++) {
for (size_t j=0; j<dirichletNodes[i].size(); j++) {
print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
}
}*/
/*for (size_t i=0; i<frictionNodes.size(); i++) {
print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
}*/
//DUNE_THROW(Dune::Exception, "Just need to stop here!");
Updaters current(
initRateUpdater(
parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunctions,
dirichletNodes,
contactNetwork.matrices(),
programState.u,
programState.v,
programState.a),
initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
programState.alpha,
nBodyAssembler.getContactCouplings(),
contactNetwork.couplings())
);
auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
/*
// set patch preconditioner for linear correction in TNNMG method
using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
PatchSolver patchSolver(gsStep,
preconditionerParset.get<size_t>("maximumIterations"),
preconditionerParset.get<double>("tolerance"),
nullptr,
preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
false); // absolute error
Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
preconditioner.setPatchSolver(patchSolver);
preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
*/
// set adaptive time stepper
typename ContactNetwork::ExternalForces externalForces;
contactNetwork.externalForces(externalForces);
auto&& noFriction = ZeroNonlinearity();
StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
externalForces, stateEnergyNorms);
UniformTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
timeStepper(stepBase, contactNetwork, current,
programState.relativeTime, programState.relativeTau);
size_t timeSteps = std::round(parset.get<double>("timeSteps.timeSteps"));
while (!timeStepper.reachedEnd()) {
programState.timeStep++;
//preconditioner.build();
iterationCount = timeStepper.advance();
programState.relativeTime = timeStepper.relativeTime_;
programState.relativeTau = timeStepper.relativeTau_;
current.rate_->extractDisplacement(programState.u);
current.rate_->extractVelocity(programState.v);
current.rate_->extractAcceleration(programState.a);
current.state_->extractAlpha(programState.alpha);
globalFriction.updateAlpha(programState.alpha);
/*print(programState.u, "current u:");
print(programState.v, "current v:");
print(programState.a, "current a:");
print(programState.alpha, "current alpha:");*/
contactNetwork.setDeformation(programState.u);
ioHandler->write(programState, contactNetwork, globalFriction, iterationCount, false);
if (programState.timeStep==timeSteps) {
std::cout << "limit of timeSteps reached!" << std::endl;
break; // TODO remove after debugging
}
if (terminationRequested) {
std::cerr << "Terminating prematurely" << std::endl;
break;
}
}
std::cout.rdbuf(coutbuf); //reset to standard output again
} catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
} catch (std::exception &e) {
std::cerr << "Standard exception: " << e.what() << std::endl;
} catch (...) {
}
ioHandler.reset();
}
src/foam/foam.cfg 0 → 100644
View file @ cdeb7bf9
# -*- mode:conf -*-
outPath = tresca # output written to ./output/outPath
gravity = 9.81 # [m/s^2]
[body0]
length = 6.0 # [m]
height = 1.0 # [m]
bulkModulus = 4.12e9 # [Pa] #2190
poissonRatio = 0.0 #0.3 # [1] #0.11
[body0.elastic]
density = 5e3 # [kg/m^3] #750
shearViscosity = 0.0 # [Pas]
bulkViscosity = 0.0 # [Pas]
[body0.viscoelastic]
density = 5e3 # [kg/m^3] #750
shearViscosity = 0.0 # [Pas]
bulkViscosity = 0.0 # [Pas]
[body1]
length = 5.00 # [m]
height = 1.00 # [m]
bulkModulus = 4.12e7 # [Pa]
poissonRatio = 0.3 # [1]
[body1.elastic]
density = 5e3 # [kg/m^3]
shearViscosity = 0.0 # [Pas]
bulkViscosity = 0.0 # [Pas]
[body1.viscoelastic]
density = 5e3 # [kg/m^3]
shearViscosity = 0.0 # [Pas]
bulkViscosity = 0.0 # [Pas]
[boundary.friction]
C = 0.0 # [Pa]
mu0 = 0.3 # [ ]
V0 = 1e-6 # [m/s]
L = 1e-5 # [m]
initialAlpha = 0.0 #-10.0 # [ ]
stateModel = AgeingLaw
frictionModel = Tresca #Tresca #None #Truncated #Regularised
[boundary.friction.weakening]
a = 0.010 # [ ]
b = 0.015 # [ ]
[boundary.friction.strengthening]
a = 0.010 # [ ]
b = 0.015 # [ ]
[boundary.neumann]
sigmaN = 0.0 # [Pa]
[boundary.dirichlet]
finalVelocity = 2e-3 #2e-4 # [m/s]
[io]
data.write = true
printProgress = true
restarts.first = 0
restarts.spacing= 50
restarts.write = true #true
vtk.write = true
[problem]
finalTime = 15 # [s] #1000
bodyCount = 2
[initialTime]
timeStep = 0
relativeTime = 0.0
relativeTau = 5e-4 # 1e-6
[timeSteps]
scheme = newmark
timeSteps = 2e3 # 2e3
[u0.solver]
maximumIterations = 100
verbosity = full
[a0.solver]
maximumIterations = 100
verbosity = full
[v.solver]
maximumIterations = 100
verbosity = quiet
[v.fpi]
maximumIterations = 10000
lambda = 0.5
[solver.tnnmg.preconditioner]
mode = additive
patchDepth = 1
maximumIterations = 2
verbosity = quiet
[solver.tnnmg.preconditioner.patchsolver]
maximumIterations = 100
verbosity = quiet
[solver.tnnmg.preconditioner.basesolver]
maximumIterations = 10000
verbosity = quiet
[solver.tnnmg.main]
pre = 1
multi = 5 # number of multigrid steps
post = 0
src/hdf5-writer.hh deleted 100644 → 0
View file @ 8fe950bc
#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
src/hdf5/frictionalboundary-writer.cc deleted 100644 → 0
View file @ 8fe950bc
#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"
src/hdf5/frictionalboundary-writer.hh deleted 100644 → 0
View file @ 8fe950bc
#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
src/hdf5/frictionalboundary-writer_tmpl.cc deleted 100644 → 0
View file @ 8fe950bc
#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);
src/hdf5/iteration-writer.cc deleted 100644 → 0
View file @ 8fe950bc
#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);
}
src/hdf5/iteration-writer.hh deleted 100644 → 0
View file @ 8fe950bc
#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
src/hdf5/patchinfo-writer.cc deleted 100644 → 0
View file @ 8fe950bc
#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"
src/hdf5/patchinfo-writer.hh deleted 100644 → 0
View file @ 8fe950bc
#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 "../one-body-problem-data/mygeometry.hh"
template <class LocalVector, class GridView> class GridEvaluator {
using Element = typename GridView::Grid::template Codim<0>::Entity;
public:
GridEvaluator(ConvexPolyhedron<LocalVector> const &weakPatch,
GridView const &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);
void write(ProgramState const &programState);
private:
HDF5::Group group_;
VertexBasis const &vertexBasis_;
GridEvaluator<LocalVector, GridView> const gridEvaluator_;
HDF5::SequenceIO<3> weakPatchGridVelocityWriter_;
};
#endif
src/hdf5/patchinfo-writer_tmpl.cc deleted 100644 → 0
View file @ 8fe950bc
#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>;
src/hdf5/restart-io.cc deleted 100644 → 0
View file @ 8fe950bc
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#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),
relativeTimeWriter_(group, "relativeTime"),
relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
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);
addEntry(relativeTimeWriter_, programState.timeStep,
programState.relativeTime);
addEntry(relativeTimeIncrementWriter_, programState.timeStep,
programState.relativeTau);
}
template <class ProgramState>
void RestartIO<ProgramState>::read(size_t timeStep,
ProgramState &programState) {
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);
readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
}
#include "restart-io_tmpl.cc"