From 0e84ee3b75029b67c3b26da44713adf83ed05375 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 11 Sep 2012 16:01:59 +0200
Subject: [PATCH] Implement Newmark scheme
---
src/enum_scheme.cc | 3 +++
src/enums.hh | 3 ++-
src/one-body-sample.cc | 55 ++++++++++++++++++++++++++++++++++++++
src/timestepping.cc | 60 ++++++++++++++++++++++++++++++++++++++++++
4 files changed, 120 insertions(+), 1 deletion(-)
diff --git a/src/enum_scheme.cc b/src/enum_scheme.cc
index 188a8645..9866d4d5 100644
--- a/src/enum_scheme.cc
+++ b/src/enum_scheme.cc
@@ -8,6 +8,9 @@ template <> struct StringToEnum<Config::scheme> {
if (s == "implicitEuler")
return Config::ImplicitEuler;
+ if (s == "newmark")
+ return Config::Newmark;
+
DUNE_THROW(Dune::Exception, "failed to parse enum");
}
};
diff --git a/src/enums.hh b/src/enums.hh
index 8e59946e..21fb6cdf 100644
--- a/src/enums.hh
+++ b/src/enums.hh
@@ -12,7 +12,8 @@ struct Config {
};
enum scheme {
ImplicitTwoStep,
- ImplicitEuler
+ ImplicitEuler,
+ Newmark
};
};
#endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 7c15b244..55b2c72e 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -47,6 +47,7 @@
#include <dune/istl/bvector.hh>
#include <dune/fufem/assemblers/functionalassembler.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>
@@ -165,6 +166,40 @@ int main(int argc, char *argv[]) {
P0Basis const p0Basis(leafView);
P1Basis const p1Basis(leafView);
+ MassAssembler<GridType, P1Basis::LocalFiniteElement,
+ P1Basis::LocalFiniteElement> const localMass;
+ MatrixType massMatrix;
+ {
+ timer.reset();
+ OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
+ .assemble(localMass, massMatrix);
+
+ // We have volume*gravity*density/area = normalstress (V*g*rho/A =
+ // sigma_n)
+ double volume = 1.0;
+ for (int i = 0; i < dim; ++i)
+ volume *= (upperRight[i] - lowerLeft[i]);
+
+ double area = 1.0;
+ for (int i = 0; i < dim; ++i)
+ if (i != 1)
+ area *= (upperRight[i] - lowerLeft[i]);
+
+ double const gravity = 9.81;
+
+ double const normalStress =
+ parset.get<double>("boundary.friction.normalstress");
+
+ // rho = sigma * A / V / g
+ // kg/m^d = N/m^(d-1) * m^(d-1) / m^d / (N/kg)
+ double const density = normalStress * area / volume / gravity;
+
+ massMatrix *= density;
+ if (parset.get<bool>("enable_timer"))
+ std::cout << "Assembled mass matrix in " << timer.elapsed() << "s"
+ << std::endl;
+ }
+
// Assemble elastic force on the body
StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const
@@ -200,6 +235,13 @@ int main(int argc, char *argv[]) {
VectorType ud(finestSize);
ud = 0.0;
+ VectorType ud_old(finestSize);
+ ud_old = 0.0;
+
+ VectorType udd(finestSize);
+ udd = 0.0;
+ VectorType udd_old(finestSize);
+ udd_old = 0.0;
SingletonVectorType alpha_old(finestSize);
alpha_old = parset.get<double>("boundary.friction.state.initial");
@@ -280,6 +322,12 @@ int main(int argc, char *argv[]) {
ell, stiffnessMatrix, u_old, u_old_old_ptr, ignoreNodes,
dirichletFunction, time, tau);
break;
+ case Config::Newmark:
+ timeSteppingScheme = Dune::make_shared<Newmark<
+ VectorType, MatrixType, decltype(dirichletFunction), dim>>(
+ ell, stiffnessMatrix, massMatrix, u_old, ud_old, udd_old,
+ ignoreNodes, dirichletFunction, time, tau);
+ break;
}
}
timeSteppingScheme->setup(problem_rhs, problem_iterate, problem_A);
@@ -307,6 +355,11 @@ int main(int argc, char *argv[]) {
timeSteppingScheme->extractSolution(problem_iterate, u);
timeSteppingScheme->extractVelocity(problem_iterate, ud);
+ udd = 0;
+ Arithmetic::addProduct(udd, 2.0 / tau, ud);
+ Arithmetic::addProduct(udd, -2.0 / tau, ud_old);
+ Arithmetic::addProduct(udd, -1.0, udd_old);
+
// Update the state
for (size_t i = 0; i < frictionalNodes.size(); ++i) {
if (frictionalNodes[i][0]) {
@@ -400,6 +453,8 @@ int main(int argc, char *argv[]) {
}
u_old_old = u_old;
u_old = u;
+ ud_old = ud;
+ udd_old = udd;
alpha_old = alpha;
// Compute von Mises stress and write everything to a file
diff --git a/src/timestepping.cc b/src/timestepping.cc
index dd9ad0e0..0a3e3dd7 100644
--- a/src/timestepping.cc
+++ b/src/timestepping.cc
@@ -1,3 +1,5 @@
+#include <dune/fufem/arithmetic.hh>
+
template <class VectorType, class MatrixType, class FunctionType, int dim>
class TimeSteppingScheme {
public:
@@ -149,3 +151,61 @@ class ImplicitTwoStep
private:
VectorType const *u_old_old_ptr;
};
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class Newmark
+ : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+ Newmark(VectorType const &_ell, MatrixType const &_A, MatrixType const &_B,
+ VectorType const &_u_old, VectorType const &_ud_old,
+ VectorType const &_udd_old,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction, double _time, double _tau)
+ : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
+ _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+ B(_B),
+ ud_old(_ud_old),
+ udd_old(_udd_old) {}
+
+ void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
+ MatrixType &problem_A) const {
+ problem_rhs = this->ell;
+ /* */ B.usmv(2.0 / this->tau, ud_old, problem_rhs);
+ /* */ B.usmv(1.0, udd_old, problem_rhs);
+ this->A.usmv(-1.0, this->u_old, problem_rhs);
+ this->A.usmv(-this->tau / 2.0, ud_old, problem_rhs);
+
+ // For fixed tau, we'd only really have to do this once
+ problem_A = this->A;
+ problem_A *= this->tau / 2.0;
+ Arithmetic::addProduct(problem_A, 2.0 / this->tau, B);
+
+ // ud_old makes a good initial iterate; we could use anything, though
+ problem_iterate = ud_old;
+
+ for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
+ if (this->dirichletNodes[i].count() == dim) {
+ problem_iterate[i] = 0;
+ this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
+ } else if (this->dirichletNodes[i][1])
+ problem_iterate[i][1] = 0; // Y direction prescribed
+ }
+
+ // u1 = u0 + tau/2 (du1 + du0)
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const {
+ solution = this->u_old;
+ Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
+ Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
+ }
+
+ void virtual extractVelocity(VectorType const &problem_iterate,
+ VectorType &velocity) const {
+ velocity = problem_iterate;
+ }
+
+private:
+ MatrixType const &B;
+ VectorType const &ud_old;
+ VectorType const &udd_old;
+};
--
GitLab