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

Do away with literate programming

parent 7537e143
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
src/.gitignore
+ 0
− 4
View file @ ab755e52
timestepping.cc
timestepping.hh
timestepping_tmpl.cc
one-body-sample.cc
\ No newline at end of file
src/Makefile.am
+ 0
− 16
View file @ ab755e52
...@@ -102,19 +102,3 @@ DISTCHECK_CONFIGURE_FLAGS = \ ...@@ -102,19 +102,3 @@ DISTCHECK_CONFIGURE_FLAGS = \
CXX="$(CXX)" CC="$(CC)" CXX="$(CXX)" CC="$(CC)"
include $(top_srcdir)/am/global-rules 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++))"
src/one-body-sample.org → src/one-body-sample.cc
+ 471
− 0
View file @ ab755e52
This diff is collapsed.
src/timestepping.cc 0 → 100644
+ 336
− 0
View file @ ab755e52
#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"
src/timestepping.hh 0 → 100644
+ 132
− 0
View file @ ab755e52
#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
src/timestepping.org deleted 100644 → 0
+ 0
− 602
View file @ 7537e143
* 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
src/timestepping_tmpl.cc 0 → 100644
+ 19
− 0
View file @ ab755e52
#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>;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment