diff --git a/src/.gitignore b/src/.gitignore
index aa8174205fc7e75281254e0edd1c8b65231201a9..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -1,4 +0,0 @@
-timestepping.cc
-timestepping.hh
-timestepping_tmpl.cc
-one-body-sample.cc
\ No newline at end of file
diff --git a/src/Makefile.am b/src/Makefile.am
index a05d4c69a6dfabfb921df57e38e4cadda6103d5d..48b85a8f3c1b653fdf626690cbf28fc379f1c93e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -102,19 +102,3 @@ DISTCHECK_CONFIGURE_FLAGS = \
CXX="$(CXX)" CC="$(CC)"
include $(top_srcdir)/am/global-rules
-
-BUILT_SOURCES = timestepping.hh timestepping.cc
-# Make sure the two are not built in parallel
-$(srcdir)/timestepping.cc: $(srcdir)/timestepping.hh
-
-EMACS ?= emacs
-
-$(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++))"
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
new file mode 100644
index 0000000000000000000000000000000000000000..624e94fc32dc3638000dc8b7ae68e190528f9d7a
--- /dev/null
+++ b/src/one-body-sample.cc
@@ -0,0 +1,471 @@
+#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 "state/stateupdater.hh"
+#include "state/ruinastateupdater.hh"
+#include "state/dieterichstateupdater.hh"
+
+int const dims = DIM;
+
+template <class VectorType, class MatrixType, class FunctionType, int dims>
+Dune::shared_ptr<TimeSteppingScheme<VectorType, MatrixType, FunctionType, dims>>
+initTimeStepper(Config::scheme scheme, FunctionType const &dirichletFunction,
+ Dune::BitSetVector<dims> 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, dims>>(
+ stiffnessMatrix, u_initial, ud_initial, ignoreNodes,
+ dirichletFunction);
+ case Config::ImplicitEuler:
+ return Dune::make_shared<
+ ImplicitEuler<VectorType, MatrixType, FunctionType, dims>>(
+ stiffnessMatrix, u_initial, ud_initial, ignoreNodes,
+ dirichletFunction);
+ case Config::Newmark:
+ return Dune::make_shared<
+ Newmark<VectorType, MatrixType, FunctionType, dims>>(
+ 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 {
+ Dune::Timer timer;
+
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree(srcdir "/one-body-sample.parset",
+ parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+
+ typedef Dune::FieldVector<double, dims> SmallVector;
+ typedef Dune::FieldMatrix<double, dims, dims> 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 mu = parset.get<double>("boundary.friction.mu");
+ auto const a = parset.get<double>("boundary.friction.a");
+ auto const b = parset.get<double>("boundary.friction.b");
+ auto const V0 = parset.get<double>("boundary.friction.V0");
+ auto const L = parset.get<double>("boundary.friction.L");
+
+ // {{{ Set up grid
+ typedef Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>
+ GridType;
+
+ Dune::FieldVector<typename GridType::ctype, dims> lowerLeft(0);
+ Dune::FieldVector<typename GridType::ctype, dims> upperRight(1);
+ upperRight[0] = parset.get<size_t>("body.width");
+ upperRight[1] = parset.get<size_t>("body.height");
+
+ Dune::array<unsigned int, dims> 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);
+
+ auto const refinements = parset.get<size_t>("grid.refinements");
+ grid->globalRefine(refinements);
+ size_t const finestSize = grid->size(grid->maxLevel(), dims);
+
+ 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);
+
+ // Set up the boundary
+ size_t specialNode = finestSize;
+ Dune::BitSetVector<dims> ignoreNodes(finestSize, false);
+ Dune::BitSetVector<1> neumannNodes(finestSize, false);
+ Dune::BitSetVector<1> frictionalNodes(finestSize, false);
+ {
+ typedef typename GridView::Codim<dims>::Iterator VertexLeafIterator;
+
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ GridView, Dune::MCMGVertexLayout> const myVertexMapper(leafView);
+
+ for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
+ assert(it->geometry().corners() == 1);
+ Dune::FieldVector<double, dims> const coordinates =
+ it->geometry().corner(0);
+ size_t const id = myVertexMapper.map(*it);
+
+ // Find the center of the lower face
+ switch (dims) {
+ 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])
+ ;
+ }
+ // Make sure that specialNode was set and points to a frictional node
+ assert(specialNode != finestSize);
+ assert(frictionalNodes[specialNode][0]);
+ };
+
+ // Set up functions for time-dependent boundary conditions
+ typedef Dune::VirtualFunction<double, double> FunctionType;
+ SharedPointerMap<std::string, FunctionType> functions;
+ initPython(functions);
+ auto const &dirichletFunction = functions.get("dirichletCondition");
+ auto const &neumannFunction = functions.get("neumannCondition");
+
+ // Set up normal stress, mass matrix, and gravity functional
+ double normalStress;
+ MatrixType massMatrix;
+ VectorType gravityFunctional;
+ {
+ double const gravity = 9.81;
+ double const density = parset.get<double>("body.density");
+
+ {
+ MassAssembler<GridType, P1Basis::LocalFiniteElement,
+ P1Basis::LocalFiniteElement> const localMass;
+ OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
+ .assemble(localMass, massMatrix);
+ massMatrix *= density;
+ };
+ {
+ double volume = 1.0;
+ for (int i = 0; i < dims; ++i)
+ volume *= (upperRight[i] - lowerLeft[i]);
+
+ double area = 1.0;
+ for (int i = 0; i < dims; ++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;
+ };
+ {
+ 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);
+ };
+ }
+ SingletonVectorType surfaceNormalStress(finestSize);
+ surfaceNormalStress = normalStress;
+
+ 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);
+
+ auto const nodalIntegrals =
+ assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
+ leafView, p1Basis, frictionalNodes);
+
+ // {{{ Initial conditions
+ 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.initial_log_state");
+ // }}}
+
+ // Set up TNNMG solver
+ auto const solver_tolerance = parset.get<double>("solver.tolerance");
+ typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
+ typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
+ MySolver<dims, MatrixType, VectorType, GridType, MyBlockProblemType>
+ mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance, *grid,
+ ignoreNodes);
+
+ Solver::VerbosityMode const verbosity =
+ parset.get<bool>("verbose") ? Solver::FULL : Solver::QUIET;
+
+ 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 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 timesteps = parset.get<size_t>("timeSteps");
+ auto const tau = parset.get<double>("endOfTime") / timesteps;
+
+ auto createRHS = [&](double _time, VectorType &_ell) {
+ assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
+ leafView, p1Basis, neumannNodes, _ell, neumannFunction, _time);
+ _ell += gravityFunctional;
+ };
+
+ auto const state_fpi_max =
+ parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
+ for (size_t run = 1; run <= timesteps; ++run) {
+ VectorType u;
+ VectorType ud;
+ SingletonVectorType alpha;
+ stateUpdater->extractState(alpha);
+
+ stateUpdater->nextTimeStep();
+ timeSteppingScheme->nextTimeStep();
+
+ double const time = tau * run;
+
+ VectorType ell(finestSize);
+ createRHS(time, ell);
+
+ 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);
+
+ auto solveDisplacementProblem = [&](VectorType &_problem_iterate,
+ SingletonVectorType const &_alpha) {
+ 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"),
+ parset.get<double>("solver.tolerance"), &energyNorm, verbosity,
+ false); // absolute error
+ overallSolver.preprocess();
+ overallSolver.solve();
+ };
+
+ // Since the velocity explodes in the quasistatic case, use the
+ // displacement as a convergence criterion
+ VectorType u_saved;
+ for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
+ solveDisplacementProblem(problem_iterate, alpha);
+
+ 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)
+ DUNE_THROW(Dune::Exception, "FPI failed to converge");
+ }
+ if (parset.get<bool>("printProgress"))
+ std::cerr << std::endl;
+ ;
+
+ 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;
+ ;
+ if (parset.get<bool>("writeVTK")) {
+ SingletonVectorType vonMisesStress;
+ 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;
+ }
+}
diff --git a/src/one-body-sample.org b/src/one-body-sample.org
deleted file mode 100644
index 9372163b20c9c3c60e139c0eec7360a41b9ce4b4..0000000000000000000000000000000000000000
--- a/src/one-body-sample.org
+++ /dev/null
@@ -1,548 +0,0 @@
-* Setup
-#+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<dims; ++i)
- volume *= (upperRight[i] - lowerLeft[i]);
-
- double area = 1.0;
- for (int i=0; i<dims; ++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++
-#+name: setupBoundary
-#+begin_src c++
- {
- typedef typename GridView::Codim<dims>::Iterator VertexLeafIterator;
-
- Dune::MultipleCodimMultipleGeomTypeMapper
- <GridView, Dune::MCMGVertexLayout> const myVertexMapper(leafView);
-
- for (auto it = leafView.begin<dims>();
- it != leafView.end<dims>();
- ++it) {
- assert(it->geometry().corners() == 1);
- Dune::FieldVector<double, dims> const coordinates = it->geometry().corner(0);
- size_t const id = myVertexMapper.map(*it);
-
- // Find the center of the lower face
- switch (dims) {
- 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])
- ;
- }
- // Make sure that specialNode was set and points to a frictional node
- assert(specialNode != finestSize);
- assert(frictionalNodes[specialNode][0]);
- }
-#+end_src
-* Writers
-#+name: createWriters
-#+begin_src c++
- 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);
-#+end_src
-#+name: writeData
-#+begin_src c++
- 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;
-#+end_src
-#+name: closeWriters
-#+begin_src c++
- state_writer.close();
- displacement_writer.close();
- velocity_writer.close();
- coefficient_writer.close();
-#+end_src
-#+name: assembleAndWriteStress
-#+begin_src c++
- if (parset.get<bool>("writeVTK")) {
- SingletonVectorType vonMisesStress;
- 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());
- }
-#+end_src
-* Misc. helpers
-#+name: printBenchmarkData
-#+begin_src c++
- if (parset.get<bool>("enableTimer"))
- std::cerr << std::endl
- << "Making " << timesteps
- << " time steps took " << timer.elapsed()
- << "s" << std::endl;
-#+end_src
-#+name: printInnerProgress
-#+begin_src c++
- if (parset.get<bool>("printProgress")) {
- std::cerr << '.';
- std::cerr.flush();
- }
-#+end_src
-#+name: printOuterProgress
-#+begin_src c++
- if (parset.get<bool>("printProgress"))
- std::cerr << std::endl;
-#+end_src
-* Functions
-#+name: defineInitTimeStepper
-#+begin_src c++
- template<class VectorType, class MatrixType, class FunctionType, int dims>
- Dune::shared_ptr
- <TimeSteppingScheme
- <VectorType,MatrixType,FunctionType,dims> >
- initTimeStepper(Config::scheme scheme,
- FunctionType const &dirichletFunction,
- Dune::BitSetVector<dims> 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,dims> >
- (stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
- case Config::ImplicitEuler:
- return Dune::make_shared<ImplicitEuler<VectorType,MatrixType,
- FunctionType,dims> >
- (stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
- case Config::Newmark:
- return Dune::make_shared <Newmark <VectorType,MatrixType,FunctionType,dims> >
- (stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
- ignoreNodes, dirichletFunction);
- }
- }
-#+end_src
-#+name: defineInitStateUpdater
-#+begin_src c++
- 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);
- }
- }
-#+end_src
-#+name: defineInitPython
-#+begin_src c++
- 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);
- }
-#+end_src
-#+name: mainBody
-#+begin_src c++
- try {
- Dune::Timer timer;
-
- Dune::ParameterTree parset;
- Dune::ParameterTreeParser::readINITree
- (srcdir "/one-body-sample.parset", parset);
- Dune::ParameterTreeParser::readOptions(argc, argv, parset);
-
- typedef Dune::FieldVector<double, dims> SmallVector;
- typedef Dune::FieldMatrix<double, dims, dims> 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 mu = parset.get<double>("boundary.friction.mu");
- auto const a = parset.get<double>("boundary.friction.a");
- auto const b = parset.get<double>("boundary.friction.b");
- auto const V0 = parset.get<double>("boundary.friction.V0");
- auto const L = parset.get<double>("boundary.friction.L");
-
- // {{{ Set up grid
- typedef Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>
- GridType;
-
- Dune::FieldVector<typename GridType::ctype, dims> lowerLeft(0);
- Dune::FieldVector<typename GridType::ctype, dims> upperRight(1);
- upperRight[0] = parset.get<size_t>("body.width");
- upperRight[1] = parset.get<size_t>("body.height");
-
- Dune::array<unsigned int, dims> 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);
-
- auto const refinements = parset.get<size_t>("grid.refinements");
- grid->globalRefine(refinements);
- size_t const finestSize = grid->size(grid->maxLevel(), dims);
-
- 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);
-
- // Set up the boundary
- size_t specialNode = finestSize;
- Dune::BitSetVector<dims> ignoreNodes(finestSize, false);
- Dune::BitSetVector<1> neumannNodes(finestSize, false);
- Dune::BitSetVector<1> frictionalNodes(finestSize, false);
- <<setupBoundary>>;
-
- // Set up functions for time-dependent boundary conditions
- typedef Dune::VirtualFunction<double, double> FunctionType;
- SharedPointerMap <std::string, FunctionType> functions;
- initPython(functions);
- auto const &dirichletFunction = functions.get("dirichletCondition");
- auto const &neumannFunction = functions.get("neumannCondition");
-
- // Set up normal stress, mass matrix, and gravity functional
- double normalStress;
- MatrixType massMatrix;
- VectorType gravityFunctional;
- {
- double const gravity = 9.81;
- double const density = parset.get<double>("body.density");
-
- <<assembleMassMatrix>>;
- <<computeNormalStress>>;
- <<computeGravitationalBodyForce>>;
- }
- SingletonVectorType surfaceNormalStress(finestSize);
- surfaceNormalStress = normalStress;
-
- MatrixType stiffnessMatrix;
- <<assembleStiffnessMatrix>>;
- EnergyNorm<MatrixType, VectorType> energyNorm(stiffnessMatrix);
-
- auto const nodalIntegrals
- = assemble_frictional<GridType, GridView, SmallVector, P1Basis>
- (leafView, p1Basis, frictionalNodes);
-
- // {{{ Initial conditions
- 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.initial_log_state");
- // }}}
-
- // Set up TNNMG solver
- auto const solver_tolerance = parset.get<double>("solver.tolerance");
- typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
- typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
- MySolver<dims, MatrixType, VectorType, GridType, MyBlockProblemType>
- mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance,
- ,*grid, ignoreNodes);
-
- Solver::VerbosityMode const verbosity
- = parset.get<bool>("verbose") ? Solver::FULL : Solver::QUIET;
-
- <<createWriters>>;
-
- 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 timesteps = parset.get<size_t>("timeSteps");
- auto const tau = parset.get<double>("endOfTime") / timesteps;
-
- auto createRHS = [&](double _time, VectorType &_ell) {
- assemble_neumann<GridType, GridView, SmallVector, P1Basis>
- (leafView, p1Basis, neumannNodes, _ell, neumannFunction, _time);
- _ell += gravityFunctional;
- };
-
- auto const state_fpi_max
- = parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
- for (size_t run = 1; run <= timesteps; ++run) {
- VectorType u;
- VectorType ud;
- SingletonVectorType alpha;
- stateUpdater->extractState(alpha);
-
- stateUpdater->nextTimeStep();
- timeSteppingScheme->nextTimeStep();
-
- double const time = tau * run;
-
- VectorType ell(finestSize);
- createRHS(time, ell);
-
- 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);
-
- auto solveDisplacementProblem = [&](VectorType &_problem_iterate,
- SingletonVectorType const &_alpha) {
- 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"),
- parset.get<double>("solver.tolerance"), &energyNorm, verbosity,
- false); // absolute error
- overallSolver.preprocess();
- overallSolver.solve();
- };
-
- // Since the velocity explodes in the quasistatic case, use the
- // displacement as a convergence criterion
- VectorType u_saved;
- for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
- solveDisplacementProblem(problem_iterate, alpha);
-
- timeSteppingScheme->postProcess(problem_iterate);
- timeSteppingScheme->extractDisplacement(u);
- timeSteppingScheme->extractVelocity(ud);
-
- stateUpdater->solve(ud);
- stateUpdater->extractState(alpha);
-
- <<printInnerProgress>>;
- 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)
- DUNE_THROW(Dune::Exception, "FPI failed to converge");
- }
- <<printOuterProgress>>;
-
- <<writeData>>;
- <<assembleAndWriteStress>>;
- }
- <<printBenchmarkData>>;
- <<closeWriters>>;
-
- 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
-* Skeleton
-#+begin_src c++ :tangle one-body-sample.cc :noweb yes
- // GENERATED -- DO NOT MODIFY
-
- #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 "state/stateupdater.hh"
- #include "state/ruinastateupdater.hh"
- #include "state/dieterichstateupdater.hh"
-
- int const dims = DIM;
-
- <<defineInitTimeStepper>>
- <<defineInitStateUpdater>>
- <<defineInitPython>>
-
- int
- main(int argc, char *argv[])
- {
- <<mainBody>>
- }
-#+end_src
diff --git a/src/timestepping.cc b/src/timestepping.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9b46235d3cfd41778a90748c306c5e8cdd8c5b32
--- /dev/null
+++ b/src/timestepping.cc
@@ -0,0 +1,336 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/arithmetic.hh>
+#include "timestepping.hh"
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::ImplicitEuler(
+ MatrixType const &_A, VectorType const &_u_initial,
+ VectorType const &_ud_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ u(_u_initial),
+ ud(_ud_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+ ud_old = ud;
+ u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::setup(
+ VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+ VectorType &problem_iterate, MatrixType &problem_A) {
+ postProcessCalled = false;
+
+ tau = _tau;
+
+ problem_rhs = ell;
+ A.mmv(u_old, problem_rhs);
+
+ // For fixed tau, we'd only really have to do this once
+ problem_A = A;
+ problem_A *= tau;
+
+ // ud_old makes a good initial iterate; we could use anything, though
+ problem_iterate = ud_old;
+
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ switch (dirichletNodes[i].count()) {
+ case 0:
+ continue;
+ case dim:
+ problem_iterate[i] = 0;
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ case 1:
+ if (dirichletNodes[i][0]) {
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ }
+ if (dirichletNodes[i][1]) {
+ problem_iterate[i][1] = 0;
+ break;
+ }
+ assert(false);
+ default:
+ assert(false);
+ }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::postProcess(
+ VectorType const &problem_iterate) {
+ postProcessCalled = true;
+
+ ud = problem_iterate;
+
+ u = u_old;
+ Arithmetic::addProduct(u, tau, ud);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType,
+ dim>::extractDisplacement(VectorType &displacement) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+ VectorType &velocity) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::clone() {
+ return new ImplicitEuler<VectorType, MatrixType, FunctionType, dim>(*this);
+}
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::ImplicitTwoStep(
+ MatrixType const &_A, VectorType const &_u_initial,
+ VectorType const &_ud_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ u(_u_initial),
+ ud(_ud_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void
+ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+ u_old_old = u_old;
+ ud_old = ud;
+ u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::setup(
+ VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+ VectorType &problem_iterate, MatrixType &problem_A) {
+ postProcessCalled = false;
+
+ tau = _tau;
+
+ switch (state) {
+ // Perform an implicit Euler step since we lack information
+ case NO_SETUP:
+ state = FIRST_SETUP;
+
+ problem_rhs = ell;
+ A.mmv(u_old, problem_rhs);
+
+ problem_A = A;
+ problem_A *= tau;
+ break;
+ case FIRST_SETUP:
+ state = SECOND_SETUP;
+ // FALLTHROUGH
+ case SECOND_SETUP:
+ problem_rhs = ell;
+ A.usmv(-4.0 / 3.0, u_old, problem_rhs);
+ A.usmv(+1.0 / 3.0, u_old_old, problem_rhs);
+
+ // For fixed tau, we'd only really have to do this once
+ problem_A = A;
+ problem_A *= 2.0 / 3.0 * tau;
+ break;
+ default:
+ assert(false);
+ }
+
+ // ud_old makes a good initial iterate; we could use anything, though
+ problem_iterate = ud_old;
+
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ switch (dirichletNodes[i].count()) {
+ case 0:
+ continue;
+ case dim:
+ problem_iterate[i] = 0;
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ case 1:
+ if (dirichletNodes[i][0]) {
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ }
+ if (dirichletNodes[i][1]) {
+ problem_iterate[i][1] = 0;
+ break;
+ }
+ assert(false);
+ default:
+ assert(false);
+ }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::postProcess(
+ VectorType const &problem_iterate) {
+ postProcessCalled = true;
+
+ ud = problem_iterate;
+
+ switch (state) {
+ case FIRST_SETUP:
+ u = u_old;
+ Arithmetic::addProduct(u, tau, ud);
+ break;
+ case SECOND_SETUP:
+ u = 0.0;
+ Arithmetic::addProduct(u, tau, ud);
+ Arithmetic::addProduct(u, 2.0, u_old);
+ Arithmetic::addProduct(u, -.5, u_old_old);
+ u *= 2.0 / 3.0;
+ break;
+ default:
+ assert(false);
+ }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitTwoStep<VectorType, MatrixType, FunctionType,
+ dim>::extractDisplacement(VectorType &displacement) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitTwoStep<VectorType, MatrixType, FunctionType,
+ dim>::extractVelocity(VectorType &velocity) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::clone() {
+ return new ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>(*this);
+}
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
+ MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+ VectorType const &_ud_initial, VectorType const &_udd_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction)
+ : A(_A),
+ B(_B),
+ u(_u_initial),
+ ud(_ud_initial),
+ udd(_udd_initial),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+ udd_old = udd;
+ ud_old = ud;
+ u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
+ VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+ VectorType &problem_iterate, MatrixType &problem_A) {
+ postProcessCalled = false;
+
+ tau = _tau;
+
+ problem_rhs = ell;
+ B.usmv(2.0 / tau, ud_old, problem_rhs);
+ B.usmv(1.0, udd_old, problem_rhs);
+ A.usmv(-1.0, u_old, problem_rhs);
+ A.usmv(-tau / 2.0, ud_old, problem_rhs);
+
+ // For fixed tau, we'd only really have to do this once
+ problem_A = A;
+ problem_A *= tau / 2.0;
+ Arithmetic::addProduct(problem_A, 2.0 / tau, B);
+
+ // ud_old makes a good initial iterate; we could use anything, though
+ problem_iterate = ud_old;
+
+ for (size_t i = 0; i < dirichletNodes.size(); ++i)
+ switch (dirichletNodes[i].count()) {
+ case 0:
+ continue;
+ case dim:
+ problem_iterate[i] = 0;
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ case 1:
+ if (dirichletNodes[i][0]) {
+ dirichletFunction.evaluate(time, problem_iterate[i][0]);
+ break;
+ }
+ if (dirichletNodes[i][1]) {
+ problem_iterate[i][1] = 0;
+ break;
+ }
+ assert(false);
+ default:
+ assert(false);
+ }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::postProcess(
+ VectorType const &problem_iterate) {
+ postProcessCalled = true;
+
+ ud = problem_iterate;
+
+ u = u_old;
+ Arithmetic::addProduct(u, tau / 2.0, ud);
+ Arithmetic::addProduct(u, tau / 2.0, ud_old);
+
+ udd = 0;
+ Arithmetic::addProduct(udd, 2.0 / tau, ud);
+ Arithmetic::addProduct(udd, -2.0 / tau, ud_old);
+ Arithmetic::addProduct(udd, -1.0, udd_old);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
+ VectorType &displacement) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+ VectorType &velocity) const {
+ if (!postProcessCalled)
+ DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+ velocity = ud;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+Newmark<VectorType, MatrixType, FunctionType, dim>::clone() {
+ return new Newmark<VectorType, MatrixType, FunctionType, dim>(*this);
+}
+
+#include "timestepping_tmpl.cc"
diff --git a/src/timestepping.hh b/src/timestepping.hh
new file mode 100644
index 0000000000000000000000000000000000000000..84da56b74f7775a079b6c3bd3dd968e992817a8f
--- /dev/null
+++ b/src/timestepping.hh
@@ -0,0 +1,132 @@
+#ifndef DUNE_TECTONIC_TIMESTEPPING_HH
+#define DUNE_TECTONIC_TIMESTEPPING_HH
+#include <dune/common/bitsetvector.hh>
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class TimeSteppingScheme {
+public:
+ void virtual nextTimeStep() = 0;
+ void virtual setup(VectorType const &ell, double _tau, double time,
+ VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A) = 0;
+
+ void virtual postProcess(VectorType const &problem_iterate) = 0;
+ void virtual extractDisplacement(VectorType &displacement) const = 0;
+ void virtual extractVelocity(VectorType &velocity) const = 0;
+
+ virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ clone() = 0;
+};
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class ImplicitEuler
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ ImplicitEuler(MatrixType const &_A, VectorType const &_u_initial,
+ VectorType const &_ud_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction);
+
+ void virtual nextTimeStep();
+ void virtual setup(VectorType const &, double, double, VectorType &,
+ VectorType &, MatrixType &);
+ void virtual postProcess(VectorType const &);
+ void virtual extractDisplacement(VectorType &) const;
+ void virtual extractVelocity(VectorType &) const;
+
+ virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ clone();
+
+private:
+ MatrixType const &A;
+ VectorType u;
+ VectorType ud;
+ Dune::BitSetVector<dim> const &dirichletNodes;
+ FunctionType const &dirichletFunction;
+
+ VectorType u_old;
+ VectorType ud_old;
+
+ double tau;
+
+ bool postProcessCalled = false;
+};
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class ImplicitTwoStep
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ ImplicitTwoStep(MatrixType const &_A, VectorType const &_u_initial,
+ VectorType const &_ud_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction);
+
+ void virtual nextTimeStep();
+ void virtual setup(VectorType const &, double, double, VectorType &,
+ VectorType &, MatrixType &);
+ void virtual postProcess(VectorType const &);
+ void virtual extractDisplacement(VectorType &) const;
+ void virtual extractVelocity(VectorType &) const;
+
+ virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ clone();
+
+private:
+ MatrixType const &A;
+ VectorType u;
+ VectorType ud;
+ Dune::BitSetVector<dim> const &dirichletNodes;
+ FunctionType const &dirichletFunction;
+
+ VectorType u_old;
+ VectorType u_old_old;
+ VectorType ud_old;
+
+ double tau;
+
+ bool postProcessCalled = false;
+
+ // Handle a lack of information
+ enum state_type {
+ NO_SETUP,
+ FIRST_SETUP,
+ SECOND_SETUP
+ };
+ state_type state = NO_SETUP;
+};
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class Newmark
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ Newmark(MatrixType const &_A, MatrixType const &_B,
+ VectorType const &_u_initial, VectorType const &_ud_initial,
+ VectorType const &_udd_initial,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction);
+
+ void virtual nextTimeStep();
+ void virtual setup(VectorType const &, double, double, VectorType &,
+ VectorType &, MatrixType &);
+ void virtual postProcess(VectorType const &);
+ void virtual extractDisplacement(VectorType &) const;
+ void virtual extractVelocity(VectorType &) const;
+
+ virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
+ clone();
+
+private:
+ MatrixType const &A;
+ MatrixType const &B;
+ VectorType u;
+ VectorType ud;
+ VectorType udd;
+ Dune::BitSetVector<dim> const &dirichletNodes;
+ FunctionType const &dirichletFunction;
+
+ VectorType u_old;
+ VectorType ud_old;
+ VectorType udd_old;
+
+ double tau;
+
+ bool postProcessCalled = false;
+};
+#endif
diff --git a/src/timestepping.org b/src/timestepping.org
deleted file mode 100644
index a234d86ebfd83e909525e6d39aa31b075228550d..0000000000000000000000000000000000000000
--- a/src/timestepping.org
+++ /dev/null
@@ -1,602 +0,0 @@
-* Preamble
-#+name: includes.hh
-#+begin_src c++
- #include <dune/common/bitsetvector.hh>
-#+end_src
-#+name: preamble.cc
-#+begin_src c++
- #ifdef HAVE_CONFIG_H
- #include "config.h"
- #endif
-
- #include <dune/fufem/arithmetic.hh>
- #include "timestepping.hh"
-#+end_src
-#+name: preamble.tex
-#+begin_src latex
- \documentclass{scrartcl}
- \usepackage{amsmath}
- \begin{document}
-#+end_src
-#+name: virtual_function_declaration
-#+begin_src c++
- void virtual nextTimeStep();
- void virtual setup(VectorType const &, double, double, VectorType &,
- VectorType &, MatrixType &);
- void virtual postProcess(VectorType const &);
- void virtual extractDisplacement(VectorType &) const;
- void virtual extractVelocity(VectorType &) const;
-
- virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
- clone();
-#+end_src
-#+name: dirichletCondition
-#+begin_src c++
- for (size_t i=0; i < dirichletNodes.size(); ++i)
- switch (dirichletNodes[i].count()) {
- case 0:
- continue;
- case dim:
- problem_iterate[i] = 0;
- dirichletFunction.evaluate(time, problem_iterate[i][0]);
- break;
- case 1:
- if (dirichletNodes[i][0]) {
- dirichletFunction.evaluate(time, problem_iterate[i][0]);
- break;
- }
- if (dirichletNodes[i][1]) {
- problem_iterate[i][1] = 0;
- break;
- }
- assert(false);
- default:
- assert(false);
- }
-#+end_src
-* Abstract TimeSteppingScheme
-#+name: TimeSteppingScheme.hh
-#+begin_src c++
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- class TimeSteppingScheme
- {
- public:
- void virtual nextTimeStep() = 0;
- void virtual
- setup(VectorType const &ell,
- double _tau,
- double time,
- VectorType &problem_rhs,
- VectorType &problem_iterate,
- MatrixType &problem_A) = 0;
-
- void virtual postProcess(VectorType const &problem_iterate) = 0;
- void virtual extractDisplacement(VectorType &displacement) const = 0;
- void virtual extractVelocity(VectorType &velocity) const = 0;
-
- virtual TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
- clone() = 0;
- };
-#+end_src
-* TimeSteppingScheme: Implicit Euler
-#+name: ImplicitEuler.hh
-#+begin_src c++ noweb: yes
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- class ImplicitEuler : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
- {
- public:
- ImplicitEuler(MatrixType const &_A,
- VectorType const &_u_initial,
- VectorType const &_ud_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction);
-
- <<virtual_function_declaration>>
-
- private:
- MatrixType const &A;
- VectorType u;
- VectorType ud;
- Dune::BitSetVector<dim> const &dirichletNodes;
- FunctionType const &dirichletFunction;
-
- VectorType u_old;
- VectorType ud_old;
-
- double tau;
-
- bool postProcessCalled = false;
- };
-#+end_src
-#+name: ImplicitEuler.cc
-#+begin_src c++
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- ImplicitEuler(MatrixType const &_A,
- VectorType const &_u_initial,
- VectorType const &_ud_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction)
- : A(_A),
- u(_u_initial),
- ud(_ud_initial),
- dirichletNodes(_dirichletNodes),
- dirichletFunction(_dirichletFunction)
- {}
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- nextTimeStep()
- {
- ud_old = ud;
- u_old = u;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- setup(VectorType const &ell,
- double _tau,
- double time,
- VectorType &problem_rhs,
- VectorType &problem_iterate,
- MatrixType &problem_A)
- {
- postProcessCalled = false;
-
- tau = _tau;
-
- problem_rhs = ell;
- A.mmv(u_old, problem_rhs);
-
- // For fixed tau, we'd only really have to do this once
- problem_A = A;
- problem_A *= tau;
-
- // ud_old makes a good initial iterate; we could use anything, though
- problem_iterate = ud_old;
-
- <<dirichletCondition>>
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- postProcess(VectorType const &problem_iterate)
- {
- postProcessCalled = true;
-
- ud = problem_iterate;
-
- u = u_old;
- Arithmetic::addProduct(u, tau, ud);
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- extractDisplacement(VectorType &displacement) const
- {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
- displacement = u;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- extractVelocity(VectorType &velocity) const
- {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
- velocity = ud;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
- ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
- clone() {
- return new ImplicitEuler<VectorType, MatrixType, FunctionType, dim>(*this);
- }
-#+end_src
-* TimeSteppingScheme: Implicit Two-Step
-#+begin_src latex :tangle twostep.tex :noweb yes
- \documentclass{scrartcl}
- \usepackage{amsmath}
- \begin{document}
- We start out with
- \begin{align}
- a(u_1, v - \dot u_1) + j(v) - j(\dot u_1) \ge \ell(v-\dot u_1)
- \end{align}
- With the two-step implicit scheme
- \begin{equation*}
- \tau \dot u_1 = \frac 32 u_1 - 2 u_0 + \frac 12 u_{-1}
- \end{equation*}
- or equivalently
- \begin{equation*}
- \frac 23 \left( \tau \dot u_1 + 2 u_0 - \frac 12 u_{-1} \right) = u_1
- \end{equation*}
- we obtain
- \begin{equation*}
- \frac 23 \tau a(\dot u_1, v - \dot u_1) + j(v) - j(\dot u_1) \ge \ell(v-\dot u_1) - a\left(\frac 43 u_0 - \frac 13 u_{-1}, v - \dot u_1\right)
- \end{equation*}
- \end{document}
-#+end_src
-#+name: ImplicitTwoStep.hh
-#+begin_src c++ noweb: yes
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- class ImplicitTwoStep : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
- {
- public:
- ImplicitTwoStep(MatrixType const &_A,
- VectorType const &_u_initial,
- VectorType const &_ud_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction);
-
- <<virtual_function_declaration>>
-
- private:
- MatrixType const &A;
- VectorType u;
- VectorType ud;
- Dune::BitSetVector<dim> const &dirichletNodes;
- FunctionType const &dirichletFunction;
-
- VectorType u_old;
- VectorType u_old_old;
- VectorType ud_old;
-
- double tau;
-
- bool postProcessCalled = false;
-
- // Handle a lack of information
- enum state_type { NO_SETUP, FIRST_SETUP, SECOND_SETUP };
- state_type state = NO_SETUP;
- };
-#+end_src
-#+name: ImplicitTwoStep.cc
-#+begin_src c++
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- ImplicitTwoStep(MatrixType const &_A,
- VectorType const &_u_initial,
- VectorType const &_ud_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction)
- : A(_A),
- u(_u_initial),
- ud(_ud_initial),
- dirichletNodes(_dirichletNodes),
- dirichletFunction(_dirichletFunction)
- {}
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- nextTimeStep()
- {
- u_old_old = u_old;
- ud_old = ud;
- u_old = u;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- setup(VectorType const &ell,
- double _tau,
- double time,
- VectorType &problem_rhs,
- VectorType &problem_iterate,
- MatrixType &problem_A)
- {
- postProcessCalled = false;
-
- tau = _tau;
-
- switch (state) {
- // Perform an implicit Euler step since we lack information
- case NO_SETUP:
- state = FIRST_SETUP;
-
- problem_rhs = ell;
- A.mmv(u_old, problem_rhs);
-
- problem_A = A;
- problem_A *= tau;
- break;
- case FIRST_SETUP:
- state = SECOND_SETUP;
- // FALLTHROUGH
- case SECOND_SETUP:
- problem_rhs = ell;
- A.usmv(-4.0/3.0, u_old, problem_rhs);
- A.usmv(+1.0/3.0, u_old_old, problem_rhs);
-
- // For fixed tau, we'd only really have to do this once
- problem_A = A;
- problem_A *= 2.0/3.0 * tau;
- break;
- default:
- assert(false);
- }
-
- // ud_old makes a good initial iterate; we could use anything, though
- problem_iterate = ud_old;
-
- <<dirichletCondition>>
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- postProcess(VectorType const &problem_iterate)
- {
- postProcessCalled = true;
-
- ud = problem_iterate;
-
- switch (state) {
- case FIRST_SETUP:
- u = u_old;
- Arithmetic::addProduct(u, tau, ud);
- break;
- case SECOND_SETUP:
- u = 0.0;
- Arithmetic::addProduct(u, tau, ud);
- Arithmetic::addProduct(u, 2.0, u_old);
- Arithmetic::addProduct(u, -.5, u_old_old);
- u *= 2.0/3.0;
- break;
- default:
- assert(false);
- }
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- extractDisplacement(VectorType &displacement) const
- {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
- displacement = u;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- extractVelocity(VectorType &velocity) const
- {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
- velocity = ud;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
- ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
- clone() {
- return new ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>(*this);
- }
-#+end_src
-* TimeSteppingScheme: Newmark
-#+begin_src latex :tangle newmark.tex :noweb yes
- <<preamble.tex>>
- \noindent The Newmark scheme in its classical form with $\gamma = 1/2$
- and $\beta = 1/4$ reads
- \begin{align}
- \label{eq:1} \dot u_1
- &= \dot u_0 + \frac \tau 2 (\ddot u_0 + \ddot u_1 )\\
- \label{eq:2} u_1
- &= u_0 + \tau \dot u_0 + \frac {\tau^2}4 ( \ddot u_0 + \ddot u_1 )\text.
- \intertext{We can also write \eqref{eq:1} as}
- \label{eq:3} \ddot u_1
- &= \frac 2\tau (\dot u_1 - \dot u_0) - \ddot u_0
- \intertext{so that it yields}
- \label{eq:4} u_1
- &= u_0 + \tau \dot u_0 + \frac {\tau^2}4 \ddot u_0 + \frac {\tau^2}4 \left(
- \frac 2\tau (\dot u_1 - \dot u_0) - \ddot u_0\right)\\
- &= u_0 + \tau \dot u_0 + \frac \tau 2 (\dot u_1 - \dot u_0)\nonumber\\
- &= u_0 + \frac \tau 2 (\dot u_1 + \dot u_0)\nonumber
- \end{align}
- in conjunction with \eqref{eq:2}. If we now consider the EVI
- \begin{align*}
- b(\ddot u_1, v - \dot u_1) + a(u_1, v - \dot u_1) + j(v) - j(\dot u_1)
- &\ge \ell (v - \dot u_1)
- \intertext{with unknowns $u_1$, $\dot u_1$, and $\ddot u_1$, we first derive}
- \frac 2\tau b(\dot u_1, v - \dot u_1)
- + a(u_1, v - \dot u_1) + j(v) - j(\dot u_1)
- &\ge \ell (v - \dot u_1) + b\left(
- \frac 2\tau \dot u_0 + \ddot u_0, v - \dot u_1\right)
- \intertext{from \eqref{eq:3} and then}
- \frac 2\tau b(\dot u_1, v - \dot u_1) + \frac \tau 2 a(\dot u_1, v - \dot u_1)
- + j(v) - j(\dot u_1)
- &\ge \ell (v - \dot u_1) + b\left(
- \frac 2\tau \dot u_0 + \ddot u_0, v - \dot u_1\right)\\
- &\quad - a\left(u_0 + \frac \tau 2 \dot u_0, v - \dot u_1\right)
- \end{align*}
- from \eqref{eq:4}. The only unknown is now $\dot u_1$.
- \end{document}
-#+end_src
-#+name: Newmark.hh
-#+begin_src c++ noweb: yes
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
- {
- public:
- Newmark(MatrixType const &_A,
- MatrixType const &_B,
- VectorType const &_u_initial,
- VectorType const &_ud_initial,
- VectorType const &_udd_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction);
-
- <<virtual_function_declaration>>
-
- private:
- MatrixType const &A;
- MatrixType const &B;
- VectorType u;
- VectorType ud;
- VectorType udd;
- Dune::BitSetVector<dim> const &dirichletNodes;
- FunctionType const &dirichletFunction;
-
- VectorType u_old;
- VectorType ud_old;
- VectorType udd_old;
-
- double tau;
-
- bool postProcessCalled = false;
- };
-#+end_src
-#+name: Newmark.cc
-#+begin_src c++
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- Newmark<VectorType, MatrixType, FunctionType, dim>::
- Newmark(MatrixType const &_A,
- MatrixType const &_B,
- VectorType const &_u_initial,
- VectorType const &_ud_initial,
- VectorType const &_udd_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction)
- : A(_A),
- B(_B),
- u(_u_initial),
- ud(_ud_initial),
- udd(_udd_initial),
- dirichletNodes(_dirichletNodes),
- dirichletFunction(_dirichletFunction)
- {}
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void Newmark<VectorType, MatrixType, FunctionType, dim>::
- nextTimeStep()
- {
- udd_old = udd;
- ud_old = ud;
- u_old = u;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void Newmark<VectorType, MatrixType, FunctionType, dim>::
- setup(VectorType const &ell,
- double _tau,
- double time,
- VectorType &problem_rhs,
- VectorType &problem_iterate,
- MatrixType &problem_A)
- {
- postProcessCalled = false;
-
- tau = _tau;
-
- problem_rhs = ell;
- B.usmv( 2.0/tau, ud_old, problem_rhs);
- B.usmv( 1.0, udd_old, problem_rhs);
- A.usmv( -1.0, u_old, problem_rhs);
- A.usmv(-tau/2.0, ud_old, problem_rhs);
-
- // For fixed tau, we'd only really have to do this once
- problem_A = A;
- problem_A *= tau/2.0;
- Arithmetic::addProduct(problem_A, 2.0/tau, B);
-
- // ud_old makes a good initial iterate; we could use anything, though
- problem_iterate = ud_old;
-
- <<dirichletCondition>>
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void Newmark<VectorType, MatrixType, FunctionType, dim>::
- postProcess(VectorType const &problem_iterate)
- {
- postProcessCalled = true;
-
- ud = problem_iterate;
-
- u = u_old;
- Arithmetic::addProduct(u, tau/2.0, ud);
- Arithmetic::addProduct(u, tau/2.0, ud_old);
-
- udd = 0;
- Arithmetic::addProduct(udd, 2.0/tau, ud);
- Arithmetic::addProduct(udd, -2.0/tau, ud_old);
- Arithmetic::addProduct(udd, -1.0, udd_old);
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void Newmark<VectorType, MatrixType, FunctionType, dim>::
- extractDisplacement(VectorType &displacement) const
- {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
- displacement = u;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- void Newmark<VectorType, MatrixType, FunctionType, dim>::
- extractVelocity(VectorType &velocity) const
- {
- if (!postProcessCalled)
- DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
- velocity = ud;
- }
-
- template <class VectorType, class MatrixType, class FunctionType, int dim>
- TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> *
- Newmark<VectorType, MatrixType, FunctionType, dim>::
- clone() {
- return new Newmark<VectorType, MatrixType, FunctionType, dim>(*this);
- }
-#+end_src
-* C++ code generation
-#+begin_src c++ :tangle timestepping.hh :noweb yes
- // GENERATED -- DO NOT MODIFY
- #ifndef DUNE_TECTONIC_TIMESTEPPING_HH
- #define DUNE_TECTONIC_TIMESTEPPING_HH
- <<includes.hh>>
-
- <<TimeSteppingScheme.hh>>
- <<ImplicitEuler.hh>>
- <<ImplicitTwoStep.hh>>
- <<Newmark.hh>>
- #endif
-#+end_src
-#+begin_src c++ :tangle timestepping.cc :noweb yes
- // GENERATED -- DO NOT MODIFY
- <<preamble.cc>>
-
- <<ImplicitEuler.cc>>
- <<ImplicitTwoStep.cc>>
- <<Newmark.cc>>
-
- #include "timestepping_tmpl.cc"
-#+end_src
-
-#+begin_src c++ :tangle timestepping_tmpl.cc :noweb yes
- // GENERATED -- DO NOT MODIFY
- #ifndef DIM
- #error DIM unset
- #endif
-
- #include <dune/common/fmatrix.hh>
- #include <dune/common/function.hh>
- #include <dune/common/fvector.hh>
- #include <dune/istl/bcrsmatrix.hh>
- #include <dune/istl/bvector.hh>
-
- 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::VirtualFunction<double, double> FunctionType;
-
- template class ImplicitEuler<VectorType, MatrixType, FunctionType, DIM>;
- template class ImplicitTwoStep<VectorType, MatrixType, FunctionType, DIM>;
- template class Newmark<VectorType, MatrixType, FunctionType, DIM>;
-#+end_src
diff --git a/src/timestepping_tmpl.cc b/src/timestepping_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3db2ba607cf744dcfedcef563ff836d4b0a2352e
--- /dev/null
+++ b/src/timestepping_tmpl.cc
@@ -0,0 +1,19 @@
+#ifndef DIM
+#error DIM unset
+#endif
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+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::VirtualFunction<double, double> FunctionType;
+
+template class ImplicitEuler<VectorType, MatrixType, FunctionType, DIM>;
+template class ImplicitTwoStep<VectorType, MatrixType, FunctionType, DIM>;
+template class Newmark<VectorType, MatrixType, FunctionType, DIM>;