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

[Output] Make all velocities/displacements relative

parent acfa70a5
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/sand-wedge.cc
+ 23
− 13
View file @ 88374d79
......@@ -234,6 +234,8 @@ int main(int argc, char *argv[]) {
u_initial = 0.0;
solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
parset.sub("u0.solver"));
Vector ur_initial(fineVertexCount);
ur_initial = 0.0;
ScalarVector alpha_initial(fineVertexCount);
alpha_initial = parset.get<double>("boundary.friction.initialLogState");
......@@ -253,6 +255,8 @@ int main(int argc, char *argv[]) {
velocityDirichletFunction.evaluate(0.0, v_initial_const);
assert(v_initial_const == 0.0);
}
Vector vr_initial(fineVertexCount);
vr_initial = 0.0;
Vector a_initial(fineVertexCount);
a_initial = 0.0;
......@@ -304,7 +308,7 @@ int main(int argc, char *argv[]) {
frictionWriter.writeKinetics(_u, _v);
frictionWriter.writeOther(c, _alpha);
};
report(u_initial, v_initial, alpha_initial);
report(ur_initial, vr_initial, alpha_initial);
MyVTKWriter<typename MyAssembler::VertexBasis,
typename MyAssembler::CellBasis> const
......@@ -314,7 +318,7 @@ int main(int argc, char *argv[]) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(
body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
}
SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities",
......@@ -334,10 +338,10 @@ int main(int argc, char *argv[]) {
std::fstream iterationWriter("iterations", std::fstream::out),
relaxationWriter("relaxation", std::fstream::out);
auto timeSteppingScheme =
initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, dirichletNodes, M, A, C,
u_initial, v_initial, a_initial);
auto timeSteppingScheme = initTimeStepper(
parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, dirichletNodes, M, A, C, u_initial,
ur_initial, v_initial, vr_initial, a_initial);
auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes,
......@@ -451,24 +455,30 @@ int main(int argc, char *argv[]) {
if (printProgress)
std::cout << std::endl;
report(u, v, alpha);
Vector ur;
Vector vr;
timeSteppingScheme->postProcessRelativeQuantities();
timeSteppingScheme->extractRelativeDisplacement(ur);
timeSteppingScheme->extractRelativeVelocity(vr);
report(ur, vr, alpha);
iterationWriter << std::endl;
relaxationWriter << std::endl;
{
BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
myAssembler.vertexBasis, v);
BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
displacement(myAssembler.vertexBasis, u);
specialVelocityWriter.write(velocity);
specialDisplacementWriter.write(displacement);
relativeVelocity(myAssembler.vertexBasis, vr);
BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
relativeDisplacement(myAssembler.vertexBasis, ur);
specialVelocityWriter.write(relativeVelocity);
specialDisplacementWriter.write(relativeDisplacement);
}
if (parset.get<bool>("io.writeVTK")) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(body.getYoungModulus(),
body.getPoissonRatio(), u, stress);
vtkWriter.write(timeStep, u, v, alpha, stress);
vtkWriter.write(timeStep, ur, vr, alpha, stress);
}
}
iterationWriter.close();
......
src/sliding-block.cc
+ 19
− 8
View file @ 88374d79
......@@ -222,6 +222,8 @@ int main(int argc, char *argv[]) {
u_initial = 0.0;
solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
parset.sub("u0.solver"));
Vector ur_initial(fineVertexCount);
ur_initial = 0.0;
ScalarVector alpha_initial(fineVertexCount);
alpha_initial =
......@@ -246,6 +248,9 @@ int main(int argc, char *argv[]) {
for (auto &x : v_initial)
x[0] = v_initial_const;
}
Vector vr_initial(fineVertexCount);
vr_initial = 0.0;
Vector a_initial(fineVertexCount);
a_initial = 0.0;
{
......@@ -285,7 +290,7 @@ int main(int argc, char *argv[]) {
frictionWriter.writeKinetics(_u, _v);
frictionWriter.writeOther(c, _alpha);
};
reportFriction(u_initial, v_initial, alpha_initial);
reportFriction(ur_initial, vr_initial, alpha_initial);
MyVTKWriter<typename MyAssembler::VertexBasis,
typename MyAssembler::CellBasis> const
......@@ -295,7 +300,7 @@ int main(int argc, char *argv[]) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(
body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
}
// Set up TNNMG solver
......@@ -310,10 +315,10 @@ int main(int argc, char *argv[]) {
std::fstream iterationWriter("iterations", std::fstream::out),
relaxationWriter("relaxation", std::fstream::out);
auto timeSteppingScheme =
initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, dirichletNodes, M, A, C,
u_initial, v_initial, a_initial);
auto timeSteppingScheme = initTimeStepper(
parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, dirichletNodes, M, A, C, u_initial,
ur_initial, v_initial, vr_initial, a_initial);
auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes,
......@@ -427,7 +432,13 @@ int main(int argc, char *argv[]) {
if (printProgress)
std::cout << std::endl;
reportFriction(u, v, alpha);
Vector ur;
Vector vr;
timeSteppingScheme->postProcessRelativeQuantities();
timeSteppingScheme->extractRelativeDisplacement(ur);
timeSteppingScheme->extractRelativeVelocity(vr);
reportFriction(ur, vr, alpha);
iterationWriter << std::endl;
relaxationWriter << std::endl;
......@@ -435,7 +446,7 @@ int main(int argc, char *argv[]) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(body.getYoungModulus(),
body.getPoissonRatio(), u, stress);
vtkWriter.write(timeStep, u, v, alpha, stress);
vtkWriter.write(timeStep, ur, vr, alpha, stress);
}
}
iterationWriter.close();
......
src/timestepping.hh
+ 10
− 5
View file @ 88374d79
......@@ -12,8 +12,11 @@ class TimeSteppingScheme {
Vector &rhs, Vector &iterate, Matrix &AB) = 0;
void virtual postProcess(Vector const &iterate) = 0;
void virtual postProcessRelativeQuantities() = 0;
void virtual extractDisplacement(Vector &displacement) const = 0;
void virtual extractRelativeDisplacement(Vector &displacement) const = 0;
void virtual extractVelocity(Vector &velocity) const = 0;
void virtual extractRelativeVelocity(Vector &velocity) const = 0;
void virtual extractOldVelocity(Vector &velocity) const = 0;
};
......@@ -27,17 +30,19 @@ initTimeStepper(Config::scheme scheme,
Dune::BitSetVector<dimension> const &velocityDirichletNodes,
Matrix const &massMatrix, Matrix const &stiffnessMatrix,
Matrix const &dampingMatrix, Vector const &u_initial,
Vector const &v_initial, Vector const &a_initial) {
Vector const &ur_initial, Vector const &v_initial,
Vector const &vr_initial, Vector const &a_initial) {
switch (scheme) {
case Config::Newmark:
return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
a_initial, velocityDirichletNodes, velocityDirichletFunction);
stiffnessMatrix, massMatrix, dampingMatrix, u_initial, ur_initial,
v_initial, vr_initial, a_initial, velocityDirichletNodes,
velocityDirichletFunction);
case Config::BackwardEuler:
return std::make_shared<
BackwardEuler<Vector, Matrix, Function, dimension>>(
stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
velocityDirichletNodes, velocityDirichletFunction);
stiffnessMatrix, massMatrix, dampingMatrix, u_initial, ur_initial,
v_initial, velocityDirichletNodes, velocityDirichletFunction);
default:
assert(false);
}
......
src/timestepping/backward_euler.cc
+ 40
− 2
View file @ 88374d79
template <class Vector, class Matrix, class Function, size_t dim>
BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
Matrix const &_A, Matrix const &_M, Matrix const &_C,
Vector const &_u_initial, Vector const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Vector const &_u_initial, Vector const &_ur_initial,
Vector const &_v_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: A(_A),
M(_M),
C(_C),
u(_u_initial),
ur(_ur_initial),
v(_v_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
......@@ -16,6 +17,7 @@ template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::nextTimeStep() {
v_o = v;
u_o = u;
ur_o = ur;
}
template <class Vector, class Matrix, class Function, size_t dim>
......@@ -23,6 +25,7 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(
Vector const &ell, double _tau, double relativeTime, Vector &rhs,
Vector &iterate, Matrix &AM) {
postProcessCalled = false;
postProcessRelativeQuantitiesCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
tau = _tau;
......@@ -87,6 +90,21 @@ void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
Arithmetic::addProduct(u, tau, v);
}
template <class Vector, class Matrix, class Function, size_t dim>
void
BackwardEuler<Vector, Matrix, Function, dim>::postProcessRelativeQuantities() {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
postProcessRelativeQuantitiesCalled = true;
vr = v;
for (size_t i = 0; i < dirichletNodes.size(); ++i)
vr[i][0] -= dirichletValue;
ur = ur_o;
Arithmetic::addProduct(ur, tau, vr);
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
......@@ -96,6 +114,16 @@ void BackwardEuler<Vector, Matrix, Function, dim>::extractDisplacement(
displacement = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractRelativeDisplacement(
Vector &displacement) const {
if (!postProcessRelativeQuantitiesCalled)
DUNE_THROW(Dune::Exception,
"It seems you forgot to call postProcessRelativeQuantities!");
displacement = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
Vector &velocity) const {
......@@ -105,6 +133,16 @@ void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
velocity = v;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::extractRelativeVelocity(
Vector &velocity) const {
if (!postProcessRelativeQuantitiesCalled)
DUNE_THROW(Dune::Exception,
"It seems you forgot to call postProcessRelativeQuantities!");
velocity = vr;
}
template <class Vector, class MatrixType, class FunctionType, size_t dim>
void BackwardEuler<Vector, MatrixType, FunctionType, dim>::extractOldVelocity(
Vector &velocity) const {
......
src/timestepping/backward_euler.hh
+ 9
− 1
View file @ 88374d79
......@@ -5,7 +5,8 @@ template <class Vector, class Matrix, class Function, size_t dim>
class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
public:
BackwardEuler(Matrix const &_A, Matrix const &_M, Matrix const &_C,
Vector const &_u_initial, Vector const &_v_initial,
Vector const &_u_initial, Vector const &_ur_initial,
Vector const &_v_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
......@@ -13,8 +14,11 @@ class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
void setup(Vector const &, double, double, Vector &, Vector &,
Matrix &) override;
void postProcess(Vector const &) override;
void postProcessRelativeQuantities() override;
void extractDisplacement(Vector &) const override;
void extractRelativeDisplacement(Vector &) const override;
void extractVelocity(Vector &) const override;
void extractRelativeVelocity(Vector &) const override;
void extractOldVelocity(Vector &) const override;
private:
......@@ -22,16 +26,20 @@ class BackwardEuler : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
Matrix const &M;
Matrix const &C;
Vector u;
Vector ur;
Vector v;
Vector vr;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o;
Vector ur_o;
Vector v_o;
double tau;
bool postProcessCalled = false;
bool postProcessRelativeQuantitiesCalled = false;
};
#endif
src/timestepping/newmark.cc
+ 43
− 1
View file @ 88374d79
template <class Vector, class Matrix, class Function, size_t dim>
Newmark<Vector, Matrix, Function, dim>::Newmark(
Matrix const &_A, Matrix const &_M, Matrix const &_C,
Vector const &_u_initial, Vector const &_v_initial,
Vector const &_u_initial, Vector const &_ur_initial,
Vector const &_v_initial, Vector const &_vr_initial,
Vector const &_a_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction)
: A(_A),
M(_M),
C(_C),
u(_u_initial),
ur(_ur_initial),
v(_v_initial),
vr(_vr_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {}
......@@ -17,7 +20,9 @@ template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::nextTimeStep() {
a_o = a;
v_o = v;
vr_o = vr;
u_o = u;
ur_o = ur;
}
template <class Vector, class Matrix, class Function, size_t dim>
......@@ -27,6 +32,7 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
Vector &rhs, Vector &iterate,
Matrix &AM) {
postProcessCalled = false;
postProcessRelativeQuantitiesCalled = false;
dirichletFunction.evaluate(relativeTime, dirichletValue);
tau = _tau;
......@@ -102,6 +108,22 @@ void Newmark<Vector, Matrix, Function, dim>::postProcess(
Arithmetic::subtractProduct(a, 1.0, a_o);
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::postProcessRelativeQuantities() {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
postProcessRelativeQuantitiesCalled = true;
vr = v;
for (size_t i = 0; i < dirichletNodes.size(); ++i)
vr[i][0] -= dirichletValue;
// u1 = tau/2 ( v1 + v0 ) + u0
ur = ur_o;
Arithmetic::addProduct(ur, tau / 2.0, vr);
Arithmetic::addProduct(ur, tau / 2.0, vr_o);
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractDisplacement(
Vector &displacement) const {
......@@ -111,6 +133,16 @@ void Newmark<Vector, Matrix, Function, dim>::extractDisplacement(
displacement = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractRelativeDisplacement(
Vector &displacement) const {
if (!postProcessRelativeQuantitiesCalled)
DUNE_THROW(Dune::Exception,
"It seems you forgot to call postProcessRelativeQuantities!");
displacement = ur;
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractVelocity(Vector &velocity)
const {
......@@ -120,6 +152,16 @@ void Newmark<Vector, Matrix, Function, dim>::extractVelocity(Vector &velocity)
velocity = v;
}
template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::extractRelativeVelocity(
Vector &velocity) const {
if (!postProcessRelativeQuantitiesCalled)
DUNE_THROW(Dune::Exception,
"It seems you forgot to call postProcessRelativeQuantities!");
velocity = vr;
}
template <class Vector, class MatrixType, class FunctionType, size_t dim>
void Newmark<Vector, MatrixType, FunctionType, dim>::extractOldVelocity(
Vector &velocity) const {
......
src/timestepping/newmark.hh
+ 10
− 1
View file @ 88374d79
......@@ -5,7 +5,8 @@ template <class Vector, class Matrix, class Function, size_t dim>
class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
public:
Newmark(Matrix const &_A, Matrix const &_M, Matrix const &_C,
Vector const &_u_initial, Vector const &_v_initial,
Vector const &_u_initial, Vector const &_ur_initial,
Vector const &_v_initial, Vector const &_vr_initial,
Vector const &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes,
Function const &_dirichletFunction);
......@@ -14,8 +15,11 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
void setup(Vector const &, double, double, Vector &, Vector &,
Matrix &) override;
void postProcess(Vector const &) override;
void postProcessRelativeQuantities() override;
void extractDisplacement(Vector &) const override;
void extractRelativeDisplacement(Vector &) const override;
void extractVelocity(Vector &) const override;
void extractRelativeVelocity(Vector &) const override;
void extractOldVelocity(Vector &) const override;
private:
......@@ -23,18 +27,23 @@ class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
Matrix const &M;
Matrix const &C;
Vector u;
Vector ur;
Vector v;
Vector vr;
Vector a;
Dune::BitSetVector<dim> const &dirichletNodes;
Function const &dirichletFunction;
double dirichletValue;
Vector u_o;
Vector ur_o;
Vector v_o;
Vector vr_o;
Vector a_o;
double tau;
bool postProcessCalled = false;
bool postProcessRelativeQuantitiesCalled = false;
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment