Skip to content
Snippets Groups Projects
Commit 06f567ac authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Transform one-body-sample into a literate program

parent dd1c1148
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/.gitignore
+ 1
− 0
View file @ 06f567ac
timestepping.cc
timestepping.hh
timestepping_tmpl.cc
one-body-sample.cc
\ No newline at end of file
src/Makefile.am
+ 5
− 0
View file @ 06f567ac
......@@ -111,3 +111,8 @@ $(srcdir)/timestepping.hh: $(srcdir)/timestepping.org
emacs -Q --batch --eval \
"(let (vc-handled-backends) \
(org-babel-tangle-file \"$<\" nil 'c++))"
$(srcdir)/one-body-sample.cc: $(srcdir)/one-body-sample.org
emacs -Q --batch --eval \
"(let (vc-handled-backends) \
(org-babel-tangle-file \"$<\" nil 'c++))"
src/one-body-sample.cc deleted 100644 → 0
+ 0
− 484
View file @ dd1c1148
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_IPOPT
#undef HAVE_IPOPT
#endif
#ifndef srcdir
#error srcdir unset
#endif
#ifndef DIM
#error DIM unset
#endif
#ifndef HAVE_PYTHON
#error Python is required
#endif
#if !HAVE_ALUGRID
#error ALUGRID is required
#endif
#include <cmath>
#include <exception>
#include <iostream>
#include <boost/format.hpp>
#include <Python.h>
#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/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/timer.hh>
#include <dune/grid/alugrid.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/sharedpointermap.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/solver.hh> // Solver::FULL
#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/myconvexproblem.hh>
#include "assemblers.hh"
#include "mysolver.hh"
#include "vtk.hh"
#include "enums.hh"
#include "enum_parser.cc"
#include "enum_state_model.cc"
#include "enum_scheme.cc"
#include "timestepping.hh"
#include "ruinastateupdater.hh"
#include "dieterichstateupdater.hh"
int const dim = DIM;
template <class GridView, class GridCorner>
void setup_boundary(GridView const &gridView,
Dune::BitSetVector<dim> &ignoreNodes,
Dune::BitSetVector<1> &neumannNodes,
Dune::BitSetVector<1> &frictionalNodes,
GridCorner const &lowerLeft, GridCorner const &upperRight,
size_t &specialNode) {
typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const myVertexMapper(gridView);
for (auto it = gridView.template begin<dim>();
it != gridView.template end<dim>(); ++it) {
assert(it->geometry().corners() == 1);
Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
size_t const id = myVertexMapper.map(*it);
// Find the center of the lower face
switch (dim) {
case 3:
if (coordinates[2] != lowerLeft[2])
break;
// fallthrough
case 2:
if (coordinates[1] == lowerLeft[1] &&
std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
specialNode = id;
break;
default:
assert(false);
}
// lower face
if (coordinates[1] == lowerLeft[1]) {
frictionalNodes[id] = true;
ignoreNodes[id][1] = true;
}
// upper face
else if (coordinates[1] == upperRight[1])
ignoreNodes[id] = true;
// right face (except for both corners)
else if (coordinates[0] == upperRight[0])
;
// left face (except for both corners)
else if (coordinates[0] == lowerLeft[0])
;
}
}
template <class VectorType, class MatrixType, class FunctionType, int dim>
Dune::shared_ptr<TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>>
initTimeStepper(Config::scheme scheme, FunctionType const &dirichletFunction,
Dune::BitSetVector<dim> const &ignoreNodes,
MatrixType const &massMatrix, MatrixType const &stiffnessMatrix,
VectorType const &u_initial, VectorType const &ud_initial,
VectorType const &udd_initial) {
switch (scheme) {
case Config::ImplicitTwoStep:
return Dune::make_shared<
ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>>(
stiffnessMatrix, u_initial, ud_initial, ignoreNodes,
dirichletFunction);
case Config::ImplicitEuler:
return Dune::make_shared<
ImplicitEuler<VectorType, MatrixType, FunctionType, dim>>(
stiffnessMatrix, u_initial, ud_initial, ignoreNodes,
dirichletFunction);
case Config::Newmark:
return Dune::make_shared<
Newmark<VectorType, MatrixType, FunctionType, dim>>(
stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
ignoreNodes, dirichletFunction);
}
}
template <class SingletonVectorType, class VectorType>
Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType>>
initStateUpdater(Config::state_model model,
SingletonVectorType const &alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes, double L) {
switch (model) {
case Config::Dieterich:
return Dune::make_shared<
DieterichStateUpdater<SingletonVectorType, VectorType>>(
alpha_initial, frictionalNodes, L);
case Config::Ruina:
return Dune::make_shared<
RuinaStateUpdater<SingletonVectorType, VectorType>>(
alpha_initial, frictionalNodes, L);
}
}
template <class FunctionMap> void initPython(FunctionMap &functions) {
Python::start();
Python::run("import sys");
Python::run("sys.path.append('" srcdir "')");
Python::import("one-body-sample")
.get("Functions")
.toC<typename FunctionMap::Base>(functions);
}
int main(int argc, char *argv[]) {
try {
typedef SharedPointerMap<std::string, Dune::VirtualFunction<double, double>>
FunctionMap;
FunctionMap functions;
initPython(functions);
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree(srcdir "/one-body-sample.parset",
parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
Dune::Timer timer;
typedef Dune::FieldVector<double, dim> SmallVector;
typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix;
typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
typedef Dune::BlockVector<SmallVector> VectorType;
typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
auto const E = parset.get<double>("body.E");
auto const nu = parset.get<double>("body.nu");
auto const solver_tolerance = parset.get<double>("solver.tolerance");
auto const refinements = parset.get<size_t>("grid.refinements");
auto const verbose = parset.get<bool>("verbose");
Solver::VerbosityMode const verbosity =
verbose ? Solver::FULL : Solver::QUIET;
// {{{ Set up grid
typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming>
GridType;
Dune::FieldVector<typename GridType::ctype, dim> lowerLeft(0);
Dune::FieldVector<typename GridType::ctype, dim> upperRight(1);
upperRight[0] = parset.get<size_t>("body.width");
upperRight[1] = parset.get<size_t>("body.height");
Dune::array<unsigned int, dim> elements;
std::fill(elements.begin(), elements.end(), 1);
elements[0] = parset.get<size_t>("body.width");
elements[1] = parset.get<size_t>("body.height");
auto grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid(
lowerLeft, upperRight, elements);
grid->globalRefine(refinements);
size_t const finestSize = grid->size(grid->maxLevel(), dim);
typedef GridType::LeafGridView GridView;
GridView const leafView = grid->leafView();
// }}}
// Set up bases
typedef P0Basis<GridView, double> P0Basis;
typedef P1NodalBasis<GridView, double> P1Basis;
P0Basis const p0Basis(leafView);
P1Basis const p1Basis(leafView);
double normalStress;
MatrixType massMatrix;
VectorType gravityFunctional;
{
double const gravity = 9.81;
double const density = parset.get<double>("body.density");
{ // Assemble mass matrix
MassAssembler<GridType, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const localMass;
OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
.assemble(localMass, massMatrix);
massMatrix *= density;
}
{ // Compute normal stress
double volume = 1.0;
for (int i = 0; i < dim; ++i)
volume *= (upperRight[i] - lowerLeft[i]);
double area = 1.0;
for (int i = 0; i < dim; ++i)
if (i != 1)
area *= (upperRight[i] - lowerLeft[i]);
// volume * gravity * density / area = normal stress
// V * g * rho / A = sigma_n
// m^d * N/kg * kg/m^d / m^(d-1) = N/m^(d-1)
normalStress = volume * gravity * density / area;
}
{ // Compute gravitational body force
SmallVector weightedGravitationalDirection(0);
weightedGravitationalDirection[1] = -density * gravity;
ConstantFunction<SmallVector, SmallVector> const gravityFunction(
weightedGravitationalDirection);
L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler(
gravityFunction);
FunctionalAssembler<P1Basis>(p1Basis)
.assemble(gravityFunctionalAssembler, gravityFunctional, true);
}
}
// Assemble elastic force on the body
MatrixType stiffnessMatrix;
{
StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const
localStiffness(E, nu);
OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
.assemble(localStiffness, stiffnessMatrix);
}
EnergyNorm<MatrixType, VectorType> energyNorm(stiffnessMatrix);
// Set up the boundary
size_t specialNode = finestSize;
Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
Dune::BitSetVector<1> neumannNodes(finestSize, false);
Dune::BitSetVector<1> frictionalNodes(finestSize, false);
setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
lowerLeft, upperRight, specialNode);
assert(specialNode != finestSize);
assert(frictionalNodes[specialNode][0]);
auto const nodalIntegrals =
assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
leafView, p1Basis, frictionalNodes);
// {{{ Initialise vectors
VectorType u(finestSize);
u = 0.0;
VectorType ud(finestSize);
ud = 0.0;
VectorType u_initial(finestSize);
u_initial = 0.0;
VectorType ud_initial(finestSize);
ud_initial = 0.0;
VectorType udd_initial(finestSize);
udd_initial = 0.0;
SingletonVectorType alpha_initial(finestSize);
alpha_initial = parset.get<double>("boundary.friction.state.initial");
SingletonVectorType alpha(alpha_initial);
SingletonVectorType vonMisesStress;
SingletonVectorType surfaceNormalStress(finestSize);
surfaceNormalStress = 0.0;
for (size_t i = 0; i < frictionalNodes.size(); ++i)
if (frictionalNodes[i][0])
surfaceNormalStress[i] = normalStress;
// }}}
typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
// Set up TNNMG solver
MySolver<dim, MatrixType, VectorType, GridType, MyBlockProblemType>
mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance, *grid,
ignoreNodes);
std::fstream state_writer("state", std::fstream::out);
std::fstream displacement_writer("displacement", std::fstream::out);
std::fstream velocity_writer("velocity", std::fstream::out);
std::fstream coefficient_writer("coefficient", std::fstream::out);
auto const L = parset.get<double>("boundary.friction.ruina.L");
auto const a = parset.get<double>("boundary.friction.ruina.a");
auto const b = parset.get<double>("boundary.friction.ruina.b");
auto const V0 = parset.get<double>("boundary.friction.V0");
auto const mu = parset.get<double>("boundary.friction.mu");
auto const timesteps = parset.get<size_t>("timeSteps");
auto const tau = parset.get<double>("endOfTime") / timesteps;
auto const &dirichletFunction = functions.get("dirichletCondition");
auto const &neumannFunction = functions.get("neumannCondition");
auto timeSteppingScheme =
initTimeStepper(parset.get<Config::scheme>("timeSteppingScheme"),
dirichletFunction, ignoreNodes, massMatrix,
stiffnessMatrix, u_initial, ud_initial, udd_initial);
auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>(
parset.get<Config::state_model>("boundary.friction.state.model"),
alpha_initial, frictionalNodes, L);
auto const state_fpi_max =
parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
for (size_t run = 1; run <= timesteps; ++run) {
double const time = tau * run;
{
stateUpdater->nextTimeStep();
VectorType ell(finestSize);
assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
ell += gravityFunctional;
// Copy everything for now
MatrixType problem_A;
VectorType problem_rhs(finestSize);
VectorType problem_iterate(finestSize);
stateUpdater->setup(tau);
timeSteppingScheme->setup(ell, tau, time, problem_rhs, problem_iterate,
problem_A);
VectorType u_saved;
for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
auto myGlobalNonlinearity =
assemble_nonlinearity<MatrixType, VectorType>(
parset.sub("boundary.friction"), *nodalIntegrals, alpha,
surfaceNormalStress);
MyConvexProblemType const myConvexProblem(
problem_A, *myGlobalNonlinearity, problem_rhs);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
auto multigridStep = mySolver.getSolver();
multigridStep->setProblem(problem_iterate, myBlockProblem);
LoopSolver<VectorType> overallSolver(
multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
solver_tolerance, &energyNorm, verbosity,
false); // absolute error
overallSolver.solve();
timeSteppingScheme->postProcess(problem_iterate);
timeSteppingScheme->extractDisplacement(u);
timeSteppingScheme->extractVelocity(ud);
stateUpdater->solve(ud);
stateUpdater->extractState(alpha);
if (parset.get<bool>("printProgress")) {
std::cerr << '.';
std::cerr.flush();
}
if (state_fpi > 1 &&
energyNorm.diff(u_saved, u) <
parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
break;
else
u_saved = u;
if (state_fpi == state_fpi_max)
std::cerr << "[ref = " << refinements
<< "]: FPI did not converge after " << state_fpi_max
<< " iterations" << std::endl;
}
if (parset.get<bool>("printProgress"))
std::cerr << std::endl;
}
{ // Write all kinds of data
state_writer << alpha[specialNode][0] << " " << std::endl;
displacement_writer << u[specialNode][0] << " " << std::endl;
velocity_writer << ud[specialNode][0] << " " << std::endl;
coefficient_writer << mu +
a *std::log(ud[specialNode].two_norm() / V0) +
b * (alpha[specialNode] + std::log(V0 / L))
<< " " << std::endl;
}
// Compute von Mises stress and write everything to a file
if (parset.get<bool>("writeVTK")) {
VonMisesStressAssembler<GridType> localStressAssembler(
E, nu,
Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>(
p1Basis, u));
FunctionalAssembler<P0Basis>(p0Basis)
.assemble(localStressAssembler, vonMisesStress, true);
writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>(
p1Basis, u, alpha, p0Basis, vonMisesStress, leafView,
(boost::format("obs%d") % run).str());
}
}
if (parset.get<bool>("enableTimer"))
std::cerr << std::endl << "Making " << timesteps << " time steps took "
<< timer.elapsed() << "s" << std::endl;
state_writer.close();
displacement_writer.close();
velocity_writer.close();
coefficient_writer.close();
Python::stop();
}
catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
}
catch (std::exception &e) {
std::cerr << "Standard exception: " << e.what() << std::endl;
}
}
src/one-body-sample.org 0 → 100644
+ 514
− 0
View file @ 06f567ac
#+name: assembleMassMatrix
#+begin_src c++
{
MassAssembler
<GridType, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const localMass;
OperatorAssembler<P1Basis,P1Basis>(p1Basis, p1Basis)
.assemble(localMass, massMatrix);
massMatrix *= density;
}
#+end_src
#+name: computeNormalStress
#+begin_src c++
{
double volume = 1.0;
for (int i=0; i<dim; ++i)
volume *= (upperRight[i] - lowerLeft[i]);
double area = 1.0;
for (int i=0; i<dim; ++i)
if (i != 1)
area *= (upperRight[i] - lowerLeft[i]);
// volume * gravity * density / area = normal stress
// V * g * rho / A = sigma_n
// m^d * N/kg * kg/m^d / m^(d-1) = N/m^(d-1)
normalStress = volume * gravity * density / area;
}
#+end_src
#+name: computeGravitationalBodyForce
#+begin_src c++
{
SmallVector weightedGravitationalDirection(0);
weightedGravitationalDirection[1] = - density * gravity;
ConstantFunction<SmallVector,SmallVector> const gravityFunction
(weightedGravitationalDirection);
L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler
(gravityFunction);
FunctionalAssembler<P1Basis>(p1Basis)
.assemble(gravityFunctionalAssembler, gravityFunctional, true);
}
#+end_src
#+name: assembleStiffnessMatrix
#+begin_src c++
{
StVenantKirchhoffAssembler
<GridType, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const localStiffness(E, nu);
OperatorAssembler<P1Basis,P1Basis>(p1Basis, p1Basis)
.assemble(localStiffness, stiffnessMatrix);
}
#+end_src c++
#+begin_src c++ :tangle one-body-sample.cc :noweb yes
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_IPOPT
#undef HAVE_IPOPT
#endif
#ifndef srcdir
#error srcdir unset
#endif
#ifndef DIM
#error DIM unset
#endif
#ifndef HAVE_PYTHON
#error Python is required
#endif
#if ! HAVE_ALUGRID
#error ALUGRID is required
#endif
#include <cmath>
#include <exception>
#include <iostream>
#include <boost/format.hpp>
#include <Python.h>
#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/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/timer.hh>
#include <dune/grid/alugrid.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/sharedpointermap.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/solver.hh> // Solver::FULL
#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/myconvexproblem.hh>
#include "assemblers.hh"
#include "mysolver.hh"
#include "vtk.hh"
#include "enums.hh"
#include "enum_parser.cc"
#include "enum_state_model.cc"
#include "enum_scheme.cc"
#include "timestepping.hh"
#include "ruinastateupdater.hh"
#include "dieterichstateupdater.hh"
int const dim = DIM;
template <class GridView, class GridCorner>
void
setup_boundary(GridView const &gridView,
Dune::BitSetVector<dim> &ignoreNodes,
Dune::BitSetVector<1> &neumannNodes,
Dune::BitSetVector<1> &frictionalNodes,
GridCorner const &lowerLeft,
GridCorner const &upperRight,
size_t &specialNode)
{
typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
Dune::MultipleCodimMultipleGeomTypeMapper
<GridView, Dune::MCMGVertexLayout> const myVertexMapper(gridView);
for (auto it = gridView.template begin<dim>();
it != gridView.template end<dim>();
++it) {
assert(it->geometry().corners() == 1);
Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
size_t const id = myVertexMapper.map(*it);
// Find the center of the lower face
switch (dim) {
case 3:
if (coordinates[2] != lowerLeft[2])
break;
// fallthrough
case 2:
if (coordinates[1] == lowerLeft[1]
&& std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
specialNode = id;
break;
default:
assert(false);
}
// lower face
if (coordinates[1] == lowerLeft[1]) {
frictionalNodes[id] = true;
ignoreNodes[id][1] = true;
}
// upper face
else if (coordinates[1] == upperRight[1])
ignoreNodes[id] = true;
// right face (except for both corners)
else if (coordinates[0] == upperRight[0])
;
// left face (except for both corners)
else if (coordinates[0] == lowerLeft[0])
;
}
}
template<class VectorType, class MatrixType, class FunctionType, int dim>
Dune::shared_ptr
<TimeSteppingScheme
<VectorType,MatrixType,FunctionType,dim> >
initTimeStepper(Config::scheme scheme,
FunctionType const &dirichletFunction,
Dune::BitSetVector<dim> const &ignoreNodes,
MatrixType const &massMatrix, MatrixType const &stiffnessMatrix,
VectorType const &u_initial,
VectorType const &ud_initial,
VectorType const &udd_initial)
{
switch (scheme) {
case Config::ImplicitTwoStep:
return Dune::make_shared<ImplicitTwoStep<VectorType,MatrixType,
FunctionType,dim> >
(stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
case Config::ImplicitEuler:
return Dune::make_shared<ImplicitEuler<VectorType,MatrixType,
FunctionType,dim> >
(stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
case Config::Newmark:
return Dune::make_shared <Newmark <VectorType,MatrixType,FunctionType,dim> >
(stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
ignoreNodes, dirichletFunction);
}
}
template<class SingletonVectorType, class VectorType>
Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType> >
initStateUpdater(Config::state_model model,
SingletonVectorType const &alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes,
double L)
{
switch (model) {
case Config::Dieterich:
return Dune::make_shared
<DieterichStateUpdater<SingletonVectorType, VectorType> >
(alpha_initial, frictionalNodes, L);
case Config::Ruina:
return Dune::make_shared
<RuinaStateUpdater<SingletonVectorType, VectorType> >
(alpha_initial, frictionalNodes, L);
}
}
template <class FunctionMap>
void
initPython(FunctionMap &functions)
{
Python::start();
Python::run("import sys");
Python::run("sys.path.append('" srcdir "')");
Python::import("one-body-sample")
.get("Functions")
.toC<typename FunctionMap::Base>(functions);
}
int
main(int argc, char *argv[])
{
try {
typedef SharedPointerMap
<std::string, Dune::VirtualFunction<double, double> > FunctionMap;
FunctionMap functions;
initPython(functions);
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree
(srcdir "/one-body-sample.parset", parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
Dune::Timer timer;
typedef Dune::FieldVector<double, dim> SmallVector;
typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix;
typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
typedef Dune::BlockVector<SmallVector> VectorType;
typedef Dune::BlockVector<Dune::FieldVector<double,1> > SingletonVectorType;
auto const E = parset.get<double>("body.E");
auto const nu = parset.get<double>("body.nu");
auto const solver_tolerance = parset.get<double>("solver.tolerance");
auto const refinements = parset.get<size_t>("grid.refinements");
auto const verbose = parset.get<bool>("verbose");
Solver::VerbosityMode const verbosity
= verbose ? Solver::FULL : Solver::QUIET;
// {{{ Set up grid
typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming>
GridType;
Dune::FieldVector<typename GridType::ctype, dim> lowerLeft(0);
Dune::FieldVector<typename GridType::ctype, dim> upperRight(1);
upperRight[0] = parset.get<size_t>("body.width");
upperRight[1] = parset.get<size_t>("body.height");
Dune::array<unsigned int, dim> elements;
std::fill(elements.begin(), elements.end(), 1);
elements[0] = parset.get<size_t>("body.width");
elements[1] = parset.get<size_t>("body.height");
auto grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid
(lowerLeft, upperRight, elements);
grid->globalRefine(refinements);
size_t const finestSize = grid->size(grid->maxLevel(), dim);
typedef GridType::LeafGridView GridView;
GridView const leafView = grid->leafView();
// }}}
// Set up bases
typedef P0Basis<GridView, double> P0Basis;
typedef P1NodalBasis<GridView, double> P1Basis;
P0Basis const p0Basis(leafView);
P1Basis const p1Basis(leafView);
double normalStress;
MatrixType massMatrix;
VectorType gravityFunctional;
{
double const gravity = 9.81;
double const density = parset.get<double>("body.density");
<<assembleMassMatrix>>;
<<computeNormalStress>>;
<<computeGravitationalBodyForce>>;
}
// Assemble elastic force on the body
MatrixType stiffnessMatrix;
<<assembleStiffnessMatrix>>;
EnergyNorm<MatrixType, VectorType> energyNorm(stiffnessMatrix);
// Set up the boundary
size_t specialNode = finestSize;
Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
Dune::BitSetVector<1> neumannNodes(finestSize, false);
Dune::BitSetVector<1> frictionalNodes(finestSize, false);
setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
lowerLeft, upperRight, specialNode);
assert(specialNode != finestSize);
assert(frictionalNodes[specialNode][0]);
auto const nodalIntegrals
= assemble_frictional<GridType, GridView, SmallVector, P1Basis>
(leafView, p1Basis, frictionalNodes);
// {{{ Initialise vectors
VectorType u(finestSize); u = 0.0;
VectorType ud(finestSize); ud = 0.0;
VectorType u_initial(finestSize); u_initial = 0.0;
VectorType ud_initial(finestSize); ud_initial = 0.0;
VectorType udd_initial(finestSize); udd_initial = 0.0;
SingletonVectorType alpha_initial(finestSize);
alpha_initial = parset.get<double>("boundary.friction.state.initial");
SingletonVectorType alpha(alpha_initial);
SingletonVectorType vonMisesStress;
SingletonVectorType surfaceNormalStress(finestSize);
surfaceNormalStress = 0.0;
for (size_t i = 0; i < frictionalNodes.size(); ++i)
if (frictionalNodes[i][0])
surfaceNormalStress[i] = normalStress;
// }}}
typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
// Set up TNNMG solver
MySolver<dim, MatrixType, VectorType, GridType, MyBlockProblemType>
mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance,
*grid, ignoreNodes);
std::fstream state_writer("state", std::fstream::out);
std::fstream displacement_writer("displacement", std::fstream::out);
std::fstream velocity_writer("velocity", std::fstream::out);
std::fstream coefficient_writer("coefficient", std::fstream::out);
auto const L = parset.get<double>("boundary.friction.ruina.L");
auto const a = parset.get<double>("boundary.friction.ruina.a");
auto const b = parset.get<double>("boundary.friction.ruina.b");
auto const V0 = parset.get<double>("boundary.friction.V0");
auto const mu = parset.get<double>("boundary.friction.mu");
auto const timesteps = parset.get<size_t>("timeSteps");
auto const tau = parset.get<double>("endOfTime") / timesteps;
auto const &dirichletFunction = functions.get("dirichletCondition");
auto const &neumannFunction = functions.get("neumannCondition");
auto timeSteppingScheme = initTimeStepper
(parset.get<Config::scheme>("timeSteppingScheme"), dirichletFunction,
ignoreNodes, massMatrix, stiffnessMatrix, u_initial, ud_initial,
udd_initial);
auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>
(parset.get<Config::state_model>("boundary.friction.state.model"),
alpha_initial, frictionalNodes, L);
auto const state_fpi_max
= parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
for (size_t run = 1; run <= timesteps; ++run) {
double const time = tau * run;
{
stateUpdater->nextTimeStep();
VectorType ell(finestSize);
assemble_neumann<GridType, GridView, SmallVector, P1Basis>
(leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
ell += gravityFunctional;
// Copy everything for now
MatrixType problem_A;
VectorType problem_rhs(finestSize);
VectorType problem_iterate(finestSize);
stateUpdater->setup(tau);
timeSteppingScheme->setup(ell, tau, time,
problem_rhs, problem_iterate, problem_A);
VectorType u_saved;
for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
auto myGlobalNonlinearity
= assemble_nonlinearity<MatrixType, VectorType>
(parset.sub("boundary.friction"), *nodalIntegrals, alpha,
surfaceNormalStress);
MyConvexProblemType const myConvexProblem
(problem_A, *myGlobalNonlinearity, problem_rhs);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
auto multigridStep = mySolver.getSolver();
multigridStep->setProblem(problem_iterate, myBlockProblem);
LoopSolver<VectorType> overallSolver
(multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
solver_tolerance, &energyNorm, verbosity, false); // absolute error
overallSolver.solve();
timeSteppingScheme->postProcess(problem_iterate);
timeSteppingScheme->extractDisplacement(u);
timeSteppingScheme->extractVelocity(ud);
stateUpdater->solve(ud);
stateUpdater->extractState(alpha);
if (parset.get<bool>("printProgress")) {
std::cerr << '.';
std::cerr.flush();
}
if (state_fpi > 1 && energyNorm.diff(u_saved, u)
< parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
break;
else
u_saved = u;
if (state_fpi == state_fpi_max)
std::cerr << "[ref = " << refinements
<< "]: FPI did not converge after "
<< state_fpi_max << " iterations" << std::endl;
}
if (parset.get<bool>("printProgress"))
std::cerr << std::endl;
}
{ // Write all kinds of data
state_writer << alpha[specialNode][0] << " " << std::endl;
displacement_writer << u[specialNode][0] << " " << std::endl;
velocity_writer << ud[specialNode][0] << " " << std::endl;
coefficient_writer << mu
+ a * std::log(ud[specialNode].two_norm() / V0)
+ b * (alpha[specialNode] + std::log(V0 / L))
<< " " << std::endl;
}
// Compute von Mises stress and write everything to a file
if (parset.get<bool>("writeVTK")) {
VonMisesStressAssembler<GridType> localStressAssembler
(E, nu,
Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>
(p1Basis, u));
FunctionalAssembler<P0Basis>(p0Basis)
.assemble(localStressAssembler, vonMisesStress, true);
writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>
(p1Basis, u, alpha, p0Basis, vonMisesStress, leafView,
(boost::format("obs%d") % run).str());
}
}
if (parset.get<bool>("enableTimer"))
std::cerr << std::endl
<< "Making " << timesteps
<< " time steps took " << timer.elapsed()
<< "s" << std::endl;
state_writer.close();
displacement_writer.close();
velocity_writer.close();
coefficient_writer.close();
Python::stop();
}
catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
}
catch (std::exception &e) {
std::cerr << "Standard exception: " << e.what() << std::endl;
}
}
#+end_src
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment