Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • podlesny/dune-tectonic
  • agnumpde/dune-tectonic
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 183 additions and 157 deletions
src/time-stepping/rate/backward_euler.hh
View file @ 09bbd1b5
...@@ -4,14 +4,15 @@ ...@@ -4,14 +4,15 @@
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> { class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
public: public:
BackwardEuler(Matrices<Matrix> const &_matrices, Vector const &_u_initial, BackwardEuler(const Matrices<Matrix> &_matrices, const std::vector<Vector> &_u_initial,
Vector const &_v_initial, Vector const &_a_initial, const std::vector<Vector> &_v_initial, const std::vector<Vector> &_a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes, const std::vector<Dune::BitSetVector<dim> > &_dirichletNodes,
Function const &_dirichletFunction); const std::vector<Function> &_dirichletFunctions);
void setup(Vector const &, double, double, Vector &, Vector &, void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
Matrix &) override; std::vector<Matrix>&) override;
void postProcess(Vector const &) override;
void postProcess(const std::vector<Vector>&) override;
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone() std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
const override; const override;
......
src/time-stepping/rate/newmark.cc
View file @ 09bbd1b5
#include <dune/solvers/common/arithmetic.hh> #include <dune/matrix-vector/axpy.hh>
#include <dune/istl/matrixindexset.hh>
#include "newmark.hh" #include "newmark.hh"
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
Newmark<Vector, Matrix, Function, dim>::Newmark( Newmark<Vector, Matrix, Function, dim>::Newmark(
Matrices<Matrix> const &_matrices, Vector const &_u_initial, const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
Vector const &_v_initial, Vector const &_a_initial, const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes, const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
Function const &_dirichletFunction) const std::vector<Function>& _dirichletFunctions)
: RateUpdater<Vector, Matrix, Function, dim>( : RateUpdater<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes, _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
_dirichletFunction) {} _dirichletFunctions) {}
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell, void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
double _tau, double _tau,
double relativeTime, double relativeTime,
Vector &rhs, Vector &iterate, std::vector<Vector>& rhs, std::vector<Vector>& iterate,
Matrix &AM) { std::vector<Matrix>& AM) {
this->dirichletFunction.evaluate(relativeTime, this->dirichletValue); for (size_t i=0; i<this->u.size(); i++) {
this->tau = _tau; this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
this->tau = _tau;
/* We start out with the formulation /* We start out with the formulation
M a + C v + A u = ell M a + C v + A u = ell
...@@ -41,58 +43,64 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell, ...@@ -41,58 +43,64 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
= [2/tau M + C + tau/2 A] v1 = [2/tau M + C + tau/2 A] v1
= ell + 2/tau M v0 + M a0 = ell + 2/tau M v0 + M a0
- tau/2 A v0 - A u0 - tau/2 A v0 - A u0
*/ */
// set up LHS (for fixed tau, we'd only really have to do this once) // set up LHS (for fixed tau, we'd only really have to do this once)
{ Matrix& LHS = AM[i];
Dune::MatrixIndexSet indices(this->matrices.elasticity.N(), {
this->matrices.elasticity.M()); Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
indices.import(this->matrices.elasticity); this->matrices.elasticity[i]->M());
indices.import(this->matrices.mass); indices.import(*this->matrices.elasticity[i]);
indices.import(this->matrices.damping); indices.import(*this->matrices.mass[i]);
indices.exportIdx(AM); indices.import(*this->matrices.damping[i]);
indices.exportIdx(LHS);
}
LHS = 0.0;
Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
// set up RHS
{
Vector& rhss = rhs[i];
rhss = ell[i];
Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
this->v_o[i]);
Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
this->v_o[i]);
Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
}
iterate = this->v_o;
const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
for (size_t k = 0; k < dirichletNodess.size(); ++k)
for (size_t j = 0; j < dim; ++j)
if (this->dirichletNodess[k][j])
iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
} }
AM = 0.0;
Arithmetic::addProduct(AM, 2.0 / this->tau, this->matrices.mass);
Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
Arithmetic::addProduct(AM, this->tau / 2.0, this->matrices.elasticity);
// set up RHS
{
rhs = ell;
Arithmetic::addProduct(rhs, 2.0 / this->tau, this->matrices.mass,
this->v_o);
Arithmetic::addProduct(rhs, this->matrices.mass, this->a_o);
Arithmetic::subtractProduct(rhs, this->tau / 2.0, this->matrices.elasticity,
this->v_o);
Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
}
iterate = this->v_o;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
for (size_t j = 0; j < dim; ++j)
if (this->dirichletNodes[i][j])
iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
} }
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void Newmark<Vector, Matrix, Function, dim>::postProcess( void Newmark<Vector, Matrix, Function, dim>::postProcess(const std::vector<Vector>& iterate) {
Vector const &iterate) {
this->postProcessCalled = true; this->postProcessCalled = true;
this->v = iterate; this->v = iterate;
// u1 = tau/2 ( v1 + v0 ) + u0 // u1 = tau/2 ( v1 + v0 ) + u0
this->u = this->u_o; this->u = this->u_o;
Arithmetic::addProduct(this->u, this->tau / 2.0, this->v);
Arithmetic::addProduct(this->u, this->tau / 2.0, this->v_o); for (size_t i=0; i<this->u.size(); i++) {
Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
// a1 = 2/tau ( v1 - v0 ) - a0 Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
this->a = 0.0;
Arithmetic::addProduct(this->a, 2.0 / this->tau, this->v); // a1 = 2/tau ( v1 - v0 ) - a0
Arithmetic::subtractProduct(this->a, 2.0 / this->tau, this->v_o); this->a[i] = 0.0;
Arithmetic::subtractProduct(this->a, 1.0, this->a_o); Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
}
} }
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
......
src/time-stepping/rate/newmark.hh
View file @ 09bbd1b5
...@@ -4,14 +4,15 @@ ...@@ -4,14 +4,15 @@
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
class Newmark : public RateUpdater<Vector, Matrix, Function, dim> { class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
public: public:
Newmark(Matrices<Matrix> const &_matrices, Vector const &_u_initial, Newmark(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
Vector const &_v_initial, Vector const &_a_initial, const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes, const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
Function const &_dirichletFunction); const std::vector<Function>& _dirichletFunctions);
void setup(Vector const &, double, double, Vector &, Vector &, void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
Matrix &) override; std::vector<Matrix>&) override;
void postProcess(Vector const &) override;
void postProcess(const std::vector<Vector>&) override;
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone() std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
const override; const override;
......
src/time-stepping/rate/rateupdater.cc
View file @ 09bbd1b5
...@@ -6,16 +6,16 @@ ...@@ -6,16 +6,16 @@
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
RateUpdater<Vector, Matrix, Function, dim>::RateUpdater( RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
Matrices<Matrix> const &_matrices, Vector const &_u_initial, const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
Vector const &_v_initial, Vector const &_a_initial, const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes, const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
Function const &_dirichletFunction) const std::vector<Function>& _dirichletFunctions)
: matrices(_matrices), : matrices(_matrices),
u(_u_initial), u(_u_initial),
v(_v_initial), v(_v_initial),
a(_a_initial), a(_a_initial),
dirichletNodes(_dirichletNodes), dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction) {} dirichletFunctions(_dirichletFunctions) {}
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() { void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
...@@ -26,17 +26,15 @@ void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() { ...@@ -26,17 +26,15 @@ void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
} }
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement( void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(std::vector<Vector>& displacements) const {
Vector &displacement) const {
if (!postProcessCalled) if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!"); DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
displacement = u; displacements = u;
} }
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity( void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(std::vector<Vector>& velocity) const {
Vector &velocity) const {
if (!postProcessCalled) if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!"); DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
...@@ -44,14 +42,12 @@ void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity( ...@@ -44,14 +42,12 @@ void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
} }
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity( void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
Vector &oldVelocity) const {
oldVelocity = v_o; oldVelocity = v_o;
} }
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration( void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(std::vector<Vector>& acceleration) const {
Vector &acceleration) const {
if (!postProcessCalled) if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!"); DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
......
src/time-stepping/rate/rateupdater.hh
View file @ 09bbd1b5
...@@ -10,10 +10,10 @@ ...@@ -10,10 +10,10 @@
template <class Vector, class Matrix, class Function, size_t dim> template <class Vector, class Matrix, class Function, size_t dim>
class RateUpdater { class RateUpdater {
protected: protected:
RateUpdater(Matrices<Matrix> const &_matrices, Vector const &_u_initial, RateUpdater(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
Vector const &_v_initial, Vector const &_a_initial, const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
Dune::BitSetVector<dim> const &_dirichletNodes, const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
Function const &_dirichletFunction); const std::vector<Function>& _dirichletFunctions);
public: public:
void nextTimeStep(); void nextTimeStep();
...@@ -21,22 +21,22 @@ class RateUpdater { ...@@ -21,22 +21,22 @@ class RateUpdater {
Vector &rhs, Vector &iterate, Matrix &AB) = 0; Vector &rhs, Vector &iterate, Matrix &AB) = 0;
void virtual postProcess(Vector const &iterate) = 0; void virtual postProcess(Vector const &iterate) = 0;
void extractDisplacement(Vector &displacement) const; void extractDisplacement(std::vector<Vector>& displacements) const;
void extractVelocity(Vector &velocity) const; void extractVelocity(std::vector<Vector>& velocity) const;
void extractOldVelocity(Vector &velocity) const; void extractOldVelocity(std::vector<Vector>& velocity) const;
void extractAcceleration(Vector &acceleration) const; void extractAcceleration(std::vector<Vector>& acceleration) const;
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone() std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone()
const = 0; const = 0;
protected: protected:
Matrices<Matrix> const &matrices; const Matrices<Matrix>& matrices;
Vector u, v, a; std::vector<Vector> u, v, a;
Dune::BitSetVector<dim> const &dirichletNodes; const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
Function const &dirichletFunction; const std::vector<Function>& dirichletFunctions;
double dirichletValue; double dirichletValue;
Vector u_o, v_o, a_o; std::vector<Vector> u_o, v_o, a_o;
double tau; double tau;
bool postProcessCalled = true; bool postProcessCalled = true;
......
src/time-stepping/rate_tmpl.cc
View file @ 09bbd1b5
...@@ -6,7 +6,10 @@ using Function = Dune::VirtualFunction<double, double>; ...@@ -6,7 +6,10 @@ using Function = Dune::VirtualFunction<double, double>;
template std::shared_ptr<RateUpdater<Vector, Matrix, Function, MY_DIM>> template std::shared_ptr<RateUpdater<Vector, Matrix, Function, MY_DIM>>
initRateUpdater<Vector, Matrix, Function, MY_DIM>( initRateUpdater<Vector, Matrix, Function, MY_DIM>(
Config::scheme scheme, Function const &velocityDirichletFunction, Config::scheme scheme,
Dune::BitSetVector<MY_DIM> const &velocityDirichletNodes, const std::vector<Function>& velocityDirichletFunctions,
Matrices<Matrix> const &matrices, Vector const &u_initial, const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
Vector const &v_initial, Vector const &a_initial); const Matrices<Matrix>& matrices,
const std::vector<Vector>& u_initial,
const std::vector<Vector>& v_initial,
const std::vector<Vector>& a_initial);
src/time-stepping/state.cc
View file @ 09bbd1b5
...@@ -7,9 +7,8 @@ ...@@ -7,9 +7,8 @@
#include "state/sliplawstateupdater.cc" #include "state/sliplawstateupdater.cc"
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater( std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
Config::stateModel model, ScalarVector const &alpha_initial, const std::vector<Dune::BitSetVector<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0) {
Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) {
switch (model) { switch (model) {
case Config::AgeingLaw: case Config::AgeingLaw:
return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>( return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
......
src/time-stepping/state.hh
View file @ 09bbd1b5
...@@ -12,7 +12,7 @@ ...@@ -12,7 +12,7 @@
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater( std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
Config::stateModel model, ScalarVector const &alpha_initial, Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes, double L, double V0); const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
#endif #endif
src/time-stepping/state/ageinglawstateupdater.cc
View file @ 09bbd1b5
...@@ -5,8 +5,10 @@ ...@@ -5,8 +5,10 @@
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater( AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes, const std::vector<ScalarVector>& _alpha_initial,
double _L, double _V0) const std::vector<Dune::BitSetVector<1>>& _nodes,
const std::vector<double>& _L,
const std::vector<double>& _V0)
: alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {} : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
...@@ -32,21 +34,28 @@ double liftSingularity(double c, double x) { ...@@ -32,21 +34,28 @@ double liftSingularity(double c, double x) {
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::solve( void AgeingLawStateUpdater<ScalarVector, Vector>::solve(
Vector const &velocity_field) { const std::vector<Vector>& velocity_field) {
for (size_t i = 0; i < nodes.size(); ++i) { for (size_t i = 0; i < alpha.size(); ++i) {
if (not toBool(nodes[i])) const auto& velocity_fieldi = velocity_field[i];
continue; const auto& nodesi = nodes[i];
auto& alphai = alpha[i];
double const V = velocity_field[i].two_norm(); auto& alpha_oi = alpha_o[i];
double const mtoL = -tau / L;
alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) + for (size_t j = 0; j < nodesi.size(); ++j) {
V0 * liftSingularity(mtoL, V)); if (not toBool(nodesi[j]))
continue;
double const V = velocity_fieldi[j].two_norm();
double const mtoL = -tau / L[i];
alphai[j] = std::log(std::exp(alpha_oi[j] + V * mtoL) +
V0[i] * liftSingularity(mtoL, V));
}
} }
} }
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha( void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
ScalarVector &_alpha) { std::vector<ScalarVector>& _alpha) {
_alpha = alpha; _alpha = alpha;
} }
......
src/time-stepping/state/ageinglawstateupdater.hh
View file @ 09bbd1b5
...@@ -6,23 +6,23 @@ ...@@ -6,23 +6,23 @@
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> { class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
public: public:
AgeingLawStateUpdater(ScalarVector const &_alpha_initial, AgeingLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
Dune::BitSetVector<1> const &_nodes, double _L, const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
double _V0); const std::vector<double>& _V0);
void nextTimeStep() override; void nextTimeStep() override;
void setup(double _tau) override; void setup(double _tau) override;
void solve(Vector const &velocity_field) override; void solve(const std::vector<Vector>& velocity_field) override;
void extractAlpha(ScalarVector &) override; void extractAlpha(std::vector<ScalarVector> &) override;
std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override; std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
private: private:
ScalarVector alpha_o; std::vector<ScalarVector> alpha_o;
ScalarVector alpha; std::vector<ScalarVector> alpha;
Dune::BitSetVector<1> const &nodes; const std::vector<Dune::BitSetVector<1>>& nodes;
double const L; const std::vector<double>& L;
double const V0; const std::vector<double>& V0;
double tau; double tau;
}; };
#endif #endif
src/time-stepping/state/sliplawstateupdater.cc
View file @ 09bbd1b5
...@@ -5,8 +5,9 @@ ...@@ -5,8 +5,9 @@
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater( SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes, const std::vector<ScalarVector>& _alpha_initial,
double _L, double _V0) const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
const std::vector<double>& _V0)
: alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {} : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
...@@ -21,21 +22,29 @@ void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) { ...@@ -21,21 +22,29 @@ void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::solve( void SlipLawStateUpdater<ScalarVector, Vector>::solve(
Vector const &velocity_field) { const std::vector<Vector>& velocity_field) {
for (size_t i = 0; i < nodes.size(); ++i) {
if (not toBool(nodes[i])) for (size_t i = 0; i < alpha.size(); ++i) {
continue; const auto& velocity_fieldi = velocity_field[i];
const auto& nodesi = nodes[i];
double const V = velocity_field[i].two_norm(); auto& alphai = alpha[i];
double const mtVoL = -tau * V / L; auto& alpha_oi = alpha_o[i];
alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) +
alpha_o[i] * std::exp(mtVoL); for (size_t j = 0; j < nodesi.size(); ++j) {
} if (not toBool(nodesi[j]))
continue;
double const V = velocity_fieldi[j].two_norm();
double const mtoL = -tau * V / L[i];
alphai[j] = (V <= 0) ? alpha_oi[j] : std::expm1(mtVoL) * std::log(V / V0[i]) +
alpha_oi[j] * std::exp(mtVoL);
}
}
} }
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha( void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
ScalarVector &_alpha) { std::vector<ScalarVector>& _alpha) {
_alpha = alpha; _alpha = alpha;
} }
......
src/time-stepping/state/sliplawstateupdater.hh
View file @ 09bbd1b5
...@@ -6,24 +6,24 @@ ...@@ -6,24 +6,24 @@
template <class ScalarVector, class Vector> template <class ScalarVector, class Vector>
class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> { class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
public: public:
SlipLawStateUpdater(ScalarVector const &_alpha_initial, SlipLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
Dune::BitSetVector<1> const &_nodes, double _L, const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
double _V0); const std::vector<double>& _V0);
void nextTimeStep() override; void nextTimeStep() override;
void setup(double _tau) override; void setup(double _tau) override;
void solve(Vector const &velocity_field) override; void solve(const std::vector<Vector>& velocity_field) override;
void extractAlpha(ScalarVector &) override; void extractAlpha(std::vector<ScalarVector> &) override;
std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override; std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
private: private:
ScalarVector alpha_o; std::vector<ScalarVector> alpha_o;
ScalarVector alpha; std::vector<ScalarVector> alpha;
Dune::BitSetVector<1> const &nodes; const std::vector<Dune::BitSetVector<1>>& nodes;
double const L; const std::vector<double>& L;
double const V0; const std::vector<double>& V0;
double tau; double tau;
}; };
#endif #endif
src/time-stepping/state/stateupdater.hh
View file @ 09bbd1b5
...@@ -7,8 +7,8 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater { ...@@ -7,8 +7,8 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
void virtual nextTimeStep() = 0; void virtual nextTimeStep() = 0;
void virtual setup(double _tau) = 0; void virtual setup(double _tau) = 0;
void virtual solve(Vector const &velocity_field) = 0; void virtual solve(const std::vector<Vector>& velocity_field) = 0;
void virtual extractAlpha(ScalarVector &alpha) = 0; void virtual extractAlpha(std::vector<ScalarVector>& alpha) = 0;
std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0; std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
}; };
......
src/time-stepping/state_tmpl.cc
View file @ 09bbd1b5
...@@ -2,5 +2,5 @@ ...@@ -2,5 +2,5 @@
template std::shared_ptr<StateUpdater<ScalarVector, Vector>> template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
initStateUpdater<ScalarVector, Vector>( initStateUpdater<ScalarVector, Vector>(
Config::stateModel model, ScalarVector const &alpha_initial, Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes, double L, double V0); const std::vector<Dune::BitSetVector<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);