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
Commits on Source (89)
Hide whitespace changes
Inline Side-by-side
Showing
with 260 additions and 717 deletions
CMakeLists.txt
View file @ 8fe950bc
...@@ -18,6 +18,8 @@ include(DuneMacros) ...@@ -18,6 +18,8 @@ include(DuneMacros)
# start a dune project with information from dune.module # start a dune project with information from dune.module
dune_project() dune_project()
find_package(HDF5 COMPONENTS C REQUIRED)
add_subdirectory("src") add_subdirectory("src")
add_subdirectory("dune") add_subdirectory("dune")
add_subdirectory("doc") add_subdirectory("doc")
......
dune-tectonic.pc.in
View file @ 8fe950bc
...@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@ ...@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
Version: @VERSION@ Version: @VERSION@
Description: dune-tectonic module Description: dune-tectonic module
URL: http://dune-project.org/ URL: http://dune-project.org/
Requires: dune-common dune-fufem dune-tnnmg Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
#Libs: -L${libdir} -ldunetectonic Libs: -L${libdir}
Cflags: -I${includedir} Cflags: -I${includedir}
dune.module
View file @ 8fe950bc
...@@ -2,9 +2,7 @@ ...@@ -2,9 +2,7 @@
# Dune module information file # # Dune module information file #
################################ ################################
#Name of the module
Module: dune-tectonic Module: dune-tectonic
Version: 2.1-1 Version: 2.5-1
Maintainer: pipping@mi.fu-berlin.de Maintainer: elias.pipping@fu-berlin.de
#depending on Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
Depends: dune-common dune-fufem dune-tnnmg
dune/tectonic/CMakeLists.txt
View file @ 8fe950bc
add_subdirectory(test)
install(FILES install(FILES
body.hh body.hh
frictiondata.hh frictiondata.hh
...@@ -12,7 +10,6 @@ install(FILES ...@@ -12,7 +10,6 @@ install(FILES
minimisation.hh minimisation.hh
myblockproblem.hh myblockproblem.hh
mydirectionalconvexfunction.hh mydirectionalconvexfunction.hh
pointtractionboundaryassembler.hh
quadraticenergy.hh quadraticenergy.hh
tectonic.hh tectonic.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
dune/tectonic/frictionpotential.hh
View file @ 8fe950bc
...@@ -29,35 +29,36 @@ class FrictionPotential { ...@@ -29,35 +29,36 @@ class FrictionPotential {
class TruncatedRateState : public FrictionPotential { class TruncatedRateState : public FrictionPotential {
public: public:
TruncatedRateState(double coefficient, double _normalStress, FrictionData _fd) TruncatedRateState(double _weight, double _weightedNormalStress,
: fd(_fd), weight(coefficient), normalStress(_normalStress) {} FrictionData _fd)
: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
double coefficientOfFriction(double V) const { double coefficientOfFriction(double V) const override {
if (V <= Vmin) if (V <= Vmin)
return 0.0; return 0.0;
return fd.a * std::log(V / Vmin); return fd.a * std::log(V / Vmin);
} }
double differential(double V) const { double differential(double V) const override {
return weight * (fd.C - normalStress * coefficientOfFriction(V)); return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
} }
double second_deriv(double V) const { double second_deriv(double V) const override {
if (V <= Vmin) if (V <= Vmin)
return 0; return 0;
return weight * (-normalStress) * (fd.a / V); return -weightedNormalStress * (fd.a / V);
} }
double regularity(double V) const { double regularity(double V) const override {
if (std::abs(V - Vmin) < 1e-14) // TODO if (std::abs(V - Vmin) < 1e-14) // TODO
return std::numeric_limits<double>::infinity(); return std::numeric_limits<double>::infinity();
return std::abs(second_deriv(V)); return std::abs(second_deriv(V));
} }
void updateAlpha(double alpha) { void updateAlpha(double alpha) override {
double const logrest = (fd.mu0 + fd.b * alpha) / fd.a; double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
Vmin = fd.V0 / std::exp(logrest); Vmin = fd.V0 / std::exp(logrest);
} }
...@@ -65,31 +66,33 @@ class TruncatedRateState : public FrictionPotential { ...@@ -65,31 +66,33 @@ class TruncatedRateState : public FrictionPotential {
private: private:
FrictionData const fd; FrictionData const fd;
double const weight; double const weight;
double const normalStress; double const weightedNormalStress;
double Vmin; double Vmin;
}; };
class RegularisedRateState : public FrictionPotential { class RegularisedRateState : public FrictionPotential {
public: public:
RegularisedRateState(double coefficient, double _normalStress, RegularisedRateState(double _weight, double _weightedNormalStress,
FrictionData _fd) FrictionData _fd)
: fd(_fd), weight(coefficient), normalStress(_normalStress) {} : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
double coefficientOfFriction(double V) const { double coefficientOfFriction(double V) const override {
return fd.a * std::asinh(0.5 * V / Vmin); return fd.a * std::asinh(0.5 * V / Vmin);
} }
double differential(double V) const { double differential(double V) const override {
return weight * (fd.C - normalStress * coefficientOfFriction(V)); return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
} }
double second_deriv(double V) const { double second_deriv(double V) const override {
return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V); return -weightedNormalStress * fd.a / std::hypot(2.0 * Vmin, V);
} }
double regularity(double V) const { return std::abs(second_deriv(V)); } double regularity(double V) const override {
return std::abs(second_deriv(V));
}
void updateAlpha(double alpha) { void updateAlpha(double alpha) override {
double const logrest = (fd.mu0 + fd.b * alpha) / fd.a; double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
Vmin = fd.V0 / std::exp(logrest); Vmin = fd.V0 / std::exp(logrest);
} }
...@@ -97,22 +100,22 @@ class RegularisedRateState : public FrictionPotential { ...@@ -97,22 +100,22 @@ class RegularisedRateState : public FrictionPotential {
private: private:
FrictionData const fd; FrictionData const fd;
double const weight; double const weight;
double const normalStress; double const weightedNormalStress;
double Vmin; double Vmin;
}; };
class TrivialFunction : public FrictionPotential { class ZeroFunction : public FrictionPotential {
public: public:
double evaluate(double) const { return 0; } double evaluate(double) const override { return 0; }
double coefficientOfFriction(double s) const { return 0; } double coefficientOfFriction(double s) const override { return 0; }
double differential(double) const { return 0; } double differential(double) const override { return 0; }
double second_deriv(double) const { return 0; } double second_deriv(double) const override { return 0; }
double regularity(double) const { return 0; } double regularity(double) const override { return 0; }
void updateAlpha(double) {} void updateAlpha(double) override {}
}; };
#endif #endif
dune/tectonic/geocoordinate.hh 0 → 100644
View file @ 8fe950bc
#ifndef DUNE_TECTONIC_GEOCOORDINATE_HH
#define DUNE_TECTONIC_GEOCOORDINATE_HH
// tiny helper to make a common piece of code pleasanter to read
template <class Geometry>
typename Geometry::GlobalCoordinate geoToPoint(Geometry geo) {
assert(geo.corners() == 1);
return geo.corner(0);
}
#endif
dune/tectonic/globalfriction.hh
View file @ 8fe950bc
...@@ -22,13 +22,12 @@ template <class Matrix, class Vector> class GlobalFriction { ...@@ -22,13 +22,12 @@ template <class Matrix, class Vector> class GlobalFriction {
using LocalMatrix = typename Matrix::block_type; using LocalMatrix = typename Matrix::block_type;
using LocalVectorType = typename Vector::block_type; using LocalVectorType = typename Vector::block_type;
size_t static const block_size = LocalVectorType::dimension; size_t static const block_size = LocalVectorType::dimension;
using Friction = LocalFriction<block_size>; using LocalNonlinearity = LocalFriction<block_size>;
double operator()(Vector const &x) const { double operator()(Vector const &x) const {
double tmp = 0; double tmp = 0;
for (size_t i = 0; i < x.size(); ++i) { for (size_t i = 0; i < x.size(); ++i) {
auto const res = restriction(i); tmp += restriction(i)(x[i]);
tmp += (*res)(x[i]);
} }
return tmp; return tmp;
} }
...@@ -36,14 +35,11 @@ template <class Matrix, class Vector> class GlobalFriction { ...@@ -36,14 +35,11 @@ template <class Matrix, class Vector> class GlobalFriction {
/* /*
Return a restriction of the outer function to the i'th node. Return a restriction of the outer function to the i'th node.
*/ */
std::shared_ptr<LocalFriction<block_size>> virtual restriction( LocalNonlinearity const virtual &restriction(size_t i) const = 0;
size_t i) const = 0;
void addHessian(Vector const &v, Matrix &hessian) const { void addHessian(Vector const &v, Matrix &hessian) const {
for (size_t i = 0; i < v.size(); ++i) { for (size_t i = 0; i < v.size(); ++i)
auto const res = restriction(i); restriction(i).addHessian(v[i], hessian[i][i]);
res->addHessian(v[i], hessian[i][i]);
}
} }
void directionalDomain(Vector const &, Vector const &, void directionalDomain(Vector const &, Vector const &,
...@@ -58,8 +54,7 @@ template <class Matrix, class Vector> class GlobalFriction { ...@@ -58,8 +54,7 @@ template <class Matrix, class Vector> class GlobalFriction {
subdifferential[0] = subdifferential[1] = 0; subdifferential[0] = subdifferential[1] = 0;
for (size_t i = 0; i < u.size(); ++i) { for (size_t i = 0; i < u.size(); ++i) {
Dune::Solvers::Interval<double> D; Dune::Solvers::Interval<double> D;
auto const res = restriction(i); restriction(i).directionalSubDiff(u[i], v[i], D);
res->directionalSubDiff(u[i], v[i], D);
subdifferential[0] += D[0]; subdifferential[0] += D[0];
subdifferential[1] += D[1]; subdifferential[1] += D[1];
} }
...@@ -71,21 +66,19 @@ template <class Matrix, class Vector> class GlobalFriction { ...@@ -71,21 +66,19 @@ template <class Matrix, class Vector> class GlobalFriction {
} }
void addGradient(Vector const &v, Vector &gradient) const { void addGradient(Vector const &v, Vector &gradient) const {
for (size_t i = 0; i < v.size(); ++i) { for (size_t i = 0; i < v.size(); ++i)
auto const res = restriction(i); restriction(i).addGradient(v[i], gradient[i]);
res->addGradient(v[i], gradient[i]);
}
} }
double regularity(size_t i, typename Vector::block_type const &x) const { double regularity(size_t i, typename Vector::block_type const &x) const {
auto const res = restriction(i); return restriction(i).regularity(x);
return res->regularity(x);
} }
void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) { ScalarVector coefficientOfFriction(Vector const &x) const {
coefficient.resize(x.size()); ScalarVector ret(x.size());
for (size_t i = 0; i < x.size(); ++i) for (size_t i = 0; i < x.size(); ++i)
coefficient[i] = restriction(i)->coefficientOfFriction(x[i]); ret[i] = restriction(i).coefficientOfFriction(x[i]);
return ret;
} }
void virtual updateAlpha(ScalarVector const &alpha) = 0; void virtual updateAlpha(ScalarVector const &alpha) = 0;
......
dune/tectonic/globalratestatefriction.hh
View file @ 8fe950bc
...@@ -10,58 +10,62 @@ ...@@ -10,58 +10,62 @@
#include <dune/istl/bcrsmatrix.hh> #include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh> #include <dune/istl/bvector.hh>
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfrictiondata.hh> #include <dune/tectonic/globalfrictiondata.hh>
#include <dune/tectonic/globalfriction.hh> #include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/index-in-sorted-range.hh>
template <class Matrix, class Vector, class ScalarFriction, class GridView> template <class Matrix, class Vector, class ScalarFriction, class GridView>
class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> { class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
public: public:
using GlobalFriction<Matrix, Vector>::block_size; using GlobalFriction<Matrix, Vector>::block_size;
using typename GlobalFriction<Matrix, Vector>::Friction; using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
private: private:
using typename GlobalFriction<Matrix, Vector>::ScalarVector; using typename GlobalFriction<Matrix, Vector>::ScalarVector;
public: public:
GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary, GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
GridView const &gridView,
GlobalFrictionData<block_size> const &frictionInfo, GlobalFrictionData<block_size> const &frictionInfo,
// Note: passing the following two makes no sense
ScalarVector const &weights, ScalarVector const &weights,
ScalarVector const &normalStress) ScalarVector const &weightedNormalStress)
: restrictions(normalStress.size()) { : restrictions(), localToGlobal(), zeroFriction() {
auto trivialNonlinearity = auto const gridView = frictionalBoundary.gridView();
std::make_shared<Friction>(std::make_shared<TrivialFunction>());
Dune::MultipleCodimMultipleGeomTypeMapper< Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView); GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
for (auto it = gridView.template begin<block_size>(); for (auto it = gridView.template begin<block_size>();
it != gridView.template end<block_size>(); ++it) { it != gridView.template end<block_size>(); ++it) {
auto const i = vertexMapper.index(*it); auto const i = vertexMapper.index(*it);
auto const coordinate = it->geometry().corner(0);
if (not frictionalBoundary.containsVertex(i)) { if (not frictionalBoundary.containsVertex(i))
restrictions[i] = trivialNonlinearity;
continue; continue;
}
auto const fp = std::make_shared<ScalarFriction>( localToGlobal.emplace_back(i);
weights[i], normalStress[i], frictionInfo(coordinate)); restrictions.emplace_back(weights[i], weightedNormalStress[i],
restrictions[i] = std::make_shared<Friction>(fp); frictionInfo(geoToPoint(it->geometry())));
} }
assert(restrictions.size() == frictionalBoundary.numVertices());
assert(localToGlobal.size() == frictionalBoundary.numVertices());
} }
void updateAlpha(ScalarVector const &alpha) override { void updateAlpha(ScalarVector const &alpha) override {
for (size_t i = 0; i < restrictions.size(); ++i) for (size_t j = 0; j < restrictions.size(); ++j)
restrictions[i]->updateAlpha(alpha[i]); restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
} }
/* /*
Return a restriction of the outer function to the i'th node. Return a restriction of the outer function to the i'th node.
*/ */
std::shared_ptr<Friction> restriction(size_t i) const override { LocalNonlinearity const &restriction(size_t i) const override {
return restrictions[i]; auto const index = indexInSortedRange(localToGlobal, i);
if (index == localToGlobal.size())
return zeroFriction;
return restrictions[index];
} }
private: private:
std::vector<std::shared_ptr<Friction>> restrictions; std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
std::vector<size_t> localToGlobal;
WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
}; };
#endif #endif
dune/tectonic/index-in-sorted-range.hh 0 → 100644
View file @ 8fe950bc
#ifndef DUNE_TECTONIC_INDEX_IN_SORTED_RANGE_HH
#define DUNE_TECTONIC_INDEX_IN_SORTED_RANGE_HH
#include <algorithm>
// returns v.size() if value does not exist
template <typename T>
size_t indexInSortedRange(std::vector<T> const &v, T value) {
size_t const specialReturnValue = v.size();
auto const b = std::begin(v);
auto const e = std::end(v);
auto const lb = std::lower_bound(b, e, value);
if (lb == e) // all elements are strictly smaller
return specialReturnValue;
if (value < *lb) // value falls between to elements
return specialReturnValue;
return std::distance(b, lb);
}
#endif
dune/tectonic/localfriction.hh
View file @ 8fe950bc
...@@ -14,39 +14,59 @@ ...@@ -14,39 +14,59 @@
template <size_t dimension> class LocalFriction { template <size_t dimension> class LocalFriction {
public: public:
virtual ~LocalFriction() {}
using VectorType = Dune::FieldVector<double, dimension>; using VectorType = Dune::FieldVector<double, dimension>;
using MatrixType = Dune::FieldMatrix<double, dimension, dimension>; using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
explicit LocalFriction(std::shared_ptr<FrictionPotential> func) void virtual updateAlpha(double alpha) = 0;
: func(func) {} double virtual regularity(VectorType const &x) const = 0;
double virtual coefficientOfFriction(VectorType const &x) const = 0;
void virtual directionalSubDiff(VectorType const &x, VectorType const &v,
Dune::Solvers::Interval<double> &D) const = 0;
double operator()(VectorType const &x) const { void virtual addHessian(VectorType const &x, MatrixType &A) const = 0;
return func->evaluate(x.two_norm());
} void virtual addGradient(VectorType const &x, VectorType &y) const = 0;
void virtual directionalDomain(
VectorType const &, VectorType const &,
Dune::Solvers::Interval<double> &dom) const = 0;
};
template <size_t dimension, class ScalarFriction>
class WrappedScalarFriction : public LocalFriction<dimension> {
using VectorType = typename LocalFriction<dimension>::VectorType;
using MatrixType = typename LocalFriction<dimension>::MatrixType;
public:
template <typename... Args>
WrappedScalarFriction(Args... args)
: func_(args...) {}
void updateAlpha(double alpha) { func->updateAlpha(alpha); } void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
double regularity(VectorType const &x) const { double regularity(VectorType const &x) const override {
double const xnorm = x.two_norm(); double const xnorm = x.two_norm();
if (xnorm <= 0.0) if (xnorm <= 0.0)
return std::numeric_limits<double>::infinity(); return std::numeric_limits<double>::infinity();
return func->regularity(xnorm); return func_.regularity(xnorm);
} }
double coefficientOfFriction(VectorType const &x) const { double coefficientOfFriction(VectorType const &x) const override {
return func->coefficientOfFriction(x.two_norm()); return func_.coefficientOfFriction(x.two_norm());
} }
// directional subdifferential: at u on the line u + t*v // directional subdifferential: at u on the line u + t*v
// u and v are assumed to be non-zero // u and v are assumed to be non-zero
void directionalSubDiff(VectorType const &x, VectorType const &v, void directionalSubDiff(VectorType const &x, VectorType const &v,
Dune::Solvers::Interval<double> &D) const { Dune::Solvers::Interval<double> &D) const override {
double const xnorm = x.two_norm(); double const xnorm = x.two_norm();
if (xnorm <= 0.0) if (xnorm <= 0.0)
D[0] = D[1] = func->differential(0.0) * v.two_norm(); D[0] = D[1] = func_.differential(0.0) * v.two_norm();
else else
D[0] = D[1] = func->differential(xnorm) * (x * v) / xnorm; D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
} }
/** Formula for the derivative: /** Formula for the derivative:
...@@ -65,14 +85,14 @@ template <size_t dimension> class LocalFriction { ...@@ -65,14 +85,14 @@ template <size_t dimension> class LocalFriction {
+ \frac {H'(|z|)}{|z|} \operatorname{id} + \frac {H'(|z|)}{|z|} \operatorname{id}
\f} \f}
*/ */
void addHessian(VectorType const &x, MatrixType &A) const { void addHessian(VectorType const &x, MatrixType &A) const override {
double const xnorm2 = x.two_norm2(); double const xnorm2 = x.two_norm2();
double const xnorm = std::sqrt(xnorm2); double const xnorm = std::sqrt(xnorm2);
if (xnorm2 <= 0.0) if (xnorm2 <= 0.0)
return; return;
double const H1 = func->differential(xnorm); double const H1 = func_.differential(xnorm);
double const H2 = func->second_deriv(xnorm); double const H2 = func_.second_deriv(xnorm);
double const tensorweight = (H2 - H1 / xnorm) / xnorm2; double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
double const idweight = H1 / xnorm; double const idweight = H1 / xnorm;
...@@ -90,22 +110,22 @@ template <size_t dimension> class LocalFriction { ...@@ -90,22 +110,22 @@ template <size_t dimension> class LocalFriction {
} }
} }
void addGradient(VectorType const &x, VectorType &y) const { void addGradient(VectorType const &x, VectorType &y) const override {
double const xnorm = x.two_norm(); double const xnorm = x.two_norm();
if (std::isinf(func->regularity(xnorm))) if (std::isinf(func_.regularity(xnorm)))
return; return;
if (xnorm > 0.0) if (xnorm > 0.0)
Arithmetic::addProduct(y, func->differential(xnorm) / xnorm, x); Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, x);
} }
void directionalDomain(VectorType const &, VectorType const &, void directionalDomain(VectorType const &, VectorType const &,
Dune::Solvers::Interval<double> &dom) const { Dune::Solvers::Interval<double> &dom) const override {
dom[0] = -std::numeric_limits<double>::max(); dom[0] = -std::numeric_limits<double>::max();
dom[1] = std::numeric_limits<double>::max(); dom[1] = std::numeric_limits<double>::max();
} }
private: private:
std::shared_ptr<FrictionPotential> const func; ScalarFriction func_;
}; };
#endif #endif
dune/tectonic/minimisation.hh
View file @ 8fe950bc
...@@ -18,7 +18,7 @@ double lineSearch(Functional const &J, ...@@ -18,7 +18,7 @@ double lineSearch(Functional const &J,
typename Functional::LocalVector const &v, typename Functional::LocalVector const &v,
Bisection const &bisection) { Bisection const &bisection) {
MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest( MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest(
J.alpha * v.two_norm2(), J.b * v, *J.phi, x, v); J.alpha * v.two_norm2(), J.b * v, J.phi, x, v);
int count; int count;
return bisection.minimize(JRest, 0.0, 0.0, count); return bisection.minimize(JRest, 0.0, 0.0, count);
} }
......
dune/tectonic/myblockproblem.hh
View file @ 8fe950bc
...@@ -4,7 +4,6 @@ ...@@ -4,7 +4,6 @@
// Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh // Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
#include <dune/common/bitsetvector.hh> #include <dune/common/bitsetvector.hh>
#include <dune/common/nullptr.hh>
#include <dune/common/parametertree.hh> #include <dune/common/parametertree.hh>
#include <dune/common/fmatrixev.hh> #include <dune/common/fmatrixev.hh>
...@@ -22,16 +21,16 @@ ...@@ -22,16 +21,16 @@
/** \brief Base class for problems where each block can be solved with a /** \brief Base class for problems where each block can be solved with a
* modified gradient method */ * modified gradient method */
template <class ConvexProblem> template <class ConvexProblem>
class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> { class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
private: private:
typedef BlockNonlinearGSProblem<ConvexProblem> BNGSP; typedef BlockNonlinearGSProblem<ConvexProblem> Base;
public: public:
using typename BNGSP::ConvexProblemType; using typename Base::ConvexProblemType;
using typename BNGSP::LocalMatrixType; using typename Base::LocalMatrixType;
using typename BNGSP::LocalVectorType; using typename Base::LocalVectorType;
using typename BNGSP::MatrixType; using typename Base::MatrixType;
using typename BNGSP::VectorType; using typename Base::VectorType;
size_t static const block_size = ConvexProblem::block_size; size_t static const block_size = ConvexProblem::block_size;
size_t static const coarse_block_size = block_size; size_t static const coarse_block_size = block_size;
...@@ -56,7 +55,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> { ...@@ -56,7 +55,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
}; };
MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem) MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
: BNGSP(parset, problem), : Base(parset, problem),
maxEigenvalues_(problem.f.size()), maxEigenvalues_(problem.f.size()),
localBisection(0.0, 1.0, 0.0, true, 0.0) { localBisection(0.0, 1.0, 0.0, true, 0.0) {
for (size_t i = 0; i < problem.f.size(); ++i) { for (size_t i = 0; i < problem.f.size(); ++i) {
...@@ -104,15 +103,10 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> { ...@@ -104,15 +103,10 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
v /= vnorm; // Rescale for numerical stability v /= vnorm; // Rescale for numerical stability
MyDirectionalConvexFunction<GlobalFriction<MatrixType, VectorType>> const auto const psi = restrict(problem_.A, problem_.f, u, v, problem_.phi);
psi(computeDirectionalA(problem_.A, v),
computeDirectionalb(problem_.A, problem_.f, u, v), problem_.phi, u,
v);
Dune::Solvers::Interval<double> D; Dune::Solvers::Interval<double> D;
psi.subDiff(0, D); psi.subDiff(0, D);
// NOTE: Numerical instability can actually get us here if (D[1] > 0) // NOTE: Numerical instability can actually get us here
if (D[1] > 0)
return 0; return 0;
int bisectionsteps = 0; int bisectionsteps = 0;
...@@ -179,17 +173,15 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> { ...@@ -179,17 +173,15 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
auto const blockEnd = std::end((*col_it)[i]); auto const blockEnd = std::end((*col_it)[i]);
for (auto blockIt = std::begin((*col_it)[i]); blockIt != blockEnd; for (auto blockIt = std::begin((*col_it)[i]); blockIt != blockEnd;
++blockIt) ++blockIt)
if (linearization.truncation[row][i] || if (linearization.truncation[row][i] or
linearization.truncation[col][blockIt.index()]) linearization.truncation[col][blockIt.index()])
*blockIt = 0.0; *blockIt = 0.0;
} }
} }
for (size_t j = 0; j < block_size; ++j) for (size_t j = 0; j < block_size; ++j)
if (linearization.truncation[row][j]) if (linearization.truncation[row][j])
linearization.b[row][j] = 0.0; linearization.b[row][j] = 0.0;
} }
for (size_t j = 0; j < block_size; ++j) for (size_t j = 0; j < block_size; ++j)
outStream << std::setw(9) << linearization.truncation.countmasked(j); outStream << std::setw(9) << linearization.truncation.countmasked(j);
} }
...@@ -203,7 +195,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> { ...@@ -203,7 +195,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
std::vector<double> maxEigenvalues_; std::vector<double> maxEigenvalues_;
// problem data // problem data
using BNGSP::problem_; using Base::problem_;
Bisection const localBisection; Bisection const localBisection;
...@@ -263,20 +255,17 @@ class MyBlockProblem<ConvexProblem>::IterateObject { ...@@ -263,20 +255,17 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
else else
ui[j] = 0; ui[j] = 0;
QuadraticEnergy<typename ConvexProblem::NonlinearityType::Friction> QuadraticEnergy<
typename ConvexProblem::NonlinearityType::LocalNonlinearity>
localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m)); localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
minimise(localJ, ui, bisection_); minimise(localJ, ui, bisection_);
} }
} }
private: private:
// problem data
ConvexProblem const &problem; ConvexProblem const &problem;
std::vector<double> maxEigenvalues_; std::vector<double> maxEigenvalues_;
Bisection const bisection_; Bisection const bisection_;
// state data for smoothing procedure used by: // state data for smoothing procedure used by:
// setIterate, updateIterate, solveLocalProblem // setIterate, updateIterate, solveLocalProblem
VectorType u; VectorType u;
......
dune/tectonic/mydirectionalconvexfunction.hh
View file @ 8fe950bc
...@@ -7,25 +7,6 @@ ...@@ -7,25 +7,6 @@
#include <dune/fufem/arithmetic.hh> #include <dune/fufem/arithmetic.hh>
#include <dune/solvers/common/interval.hh> #include <dune/solvers/common/interval.hh>
/*
1/2 <A(u + hv),u + hv> - <b, u + hv>
= 1/2 <Av,v> h^2 - <b - Au, v> h + const.
localA = <Av,v>
localb = <b - Au, v>
*/
template <class Matrix, class Vector>
double computeDirectionalA(Matrix const &A, Vector const &v) {
return Arithmetic::Axy(A, v, v);
}
template <class Matrix, class Vector>
double computeDirectionalb(Matrix const &A, Vector const &b, Vector const &u,
Vector const &v) {
return Arithmetic::bmAxy(A, b, u, v);
}
template <class Nonlinearity> class MyDirectionalConvexFunction { template <class Nonlinearity> class MyDirectionalConvexFunction {
public: public:
using Vector = typename Nonlinearity::VectorType; using Vector = typename Nonlinearity::VectorType;
...@@ -65,4 +46,23 @@ template <class Nonlinearity> class MyDirectionalConvexFunction { ...@@ -65,4 +46,23 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
Dune::Solvers::Interval<double> dom; Dune::Solvers::Interval<double> dom;
}; };
/*
1/2 <A(u + hv),u + hv> - <b, u + hv>
= 1/2 <Av,v> h^2 - <b - Au, v> h + const.
localA = <Av,v>
localb = <b - Au, v>
*/
template <class Matrix, class Vector, class Nonlinearity>
MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
Vector const &b,
Vector const &u,
Vector const &v,
Nonlinearity const &phi) {
return MyDirectionalConvexFunction<Nonlinearity>(
Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
}
#endif #endif
dune/tectonic/polyhedrondistance.hh deleted 100644 → 0
View file @ 4bea0810
#ifndef DUNE_TECTONIC_POLYHEDRONDISTANCE_HH
#define DUNE_TECTONIC_POLYHEDRONDISTANCE_HH
// Based on the closest point projection from dune-contact
#include <dune/common/dynmatrix.hh>
#include <dune/solvers/common/arithmetic.hh>
template <class Coordinate> struct ConvexPolyhedron {
std::vector<Coordinate> vertices;
template <class DynamicVector>
Coordinate barycentricToGlobal(DynamicVector const &b) const {
assert(b.size() == vertices.size());
Coordinate r(0);
for (size_t i = 0; i < b.size(); i++)
Arithmetic::addProduct(r, b[i], vertices[i]);
return r;
}
template <class DynamicVector>
void sanityCheck(DynamicVector const &b, double tol) const {
double sum = 0;
for (size_t i = 0; i < b.size(); i++) {
if (b[i] < -tol)
DUNE_THROW(Dune::Exception, "No barycentric coords " << b[1] << " nr. "
<< i);
sum += b[i];
}
if (sum > 1 + tol)
DUNE_THROW(Dune::Exception, "No barycentric coords, sum: " << sum);
}
};
template <class Coordinate, class Matrix>
void populateMatrix(ConvexPolyhedron<Coordinate> const &segment,
Matrix &matrix) {
size_t const nCorners = segment.vertices.size();
for (size_t i = 0; i < nCorners; i++)
for (size_t j = 0; j <= i; j++)
matrix[i][j] = segment.vertices[i] * segment.vertices[j];
// make symmetric
for (size_t i = 0; i < nCorners; i++)
for (size_t j = i + 1; j < nCorners; j++)
matrix[i][j] = matrix[j][i];
}
/**
* The functional to be minimized is
*
* 0.5*norm(\sum_i lambda_i corner_i - target)^2
*
* where lambda_i are barycentric coordinates, i.e.
*
* 0 \le lambda_i \le 1 and \sum_i lambda_i=1
*
* The resulting quadratic minimization problem is given by
*
* 0.5 lambda^T*A*lambda - l^T*lambda
*
* with A_ij = (corner_i,corner_j) and l_i = (corner_i, target)
*/
template <class Coordinate>
void approximate(Coordinate const &target, Dune::DynamicMatrix<double> const &A,
Dune::DynamicVector<double> const &l,
ConvexPolyhedron<Coordinate> const &segment,
Dune::DynamicVector<double> &sol) {
size_t nCorners = segment.vertices.size();
// compute the residual
Dune::DynamicVector<double> rhs = l;
// compute rhs -= A*sol_
for (size_t k = 0; k < nCorners; k++)
for (size_t j = 0; j < nCorners; j++)
rhs[k] -= A[k][j] * sol[j];
// use the edge vectors as search directions
double alpha = 0.0;
for (size_t j1 = 0; j1 < nCorners - 1; ++j1) {
for (size_t j2 = j1 + 1; j2 < nCorners; ++j2) {
// compute matrix entry and rhs for edge direction
double a = A[j1][j1] - A[j1][j2] - A[j2][j1] + A[j2][j2];
double b = rhs[j1] - rhs[j2];
// compute minimizer for correction
if (a > 0) {
alpha = b / a;
double const largestAlpha = sol[j2];
double const smallestAlpha = -sol[j1];
// project alpha such that we stay positive
if (alpha > largestAlpha)
alpha = largestAlpha;
else if (alpha < smallestAlpha)
alpha = smallestAlpha;
} else {
// if the quadratic term vanishes, we have a linear function, which
// attains its extrema on the boundary
double sum = sol[j1] + sol[j2];
double lValue = 0.5 * A[j1][j1] * sum * sum - rhs[j1] * sum;
double uValue = 0.5 * A[j2][j2] * sum * sum - rhs[j2] * sum;
alpha = lValue < uValue ? sol[j2] : -sol[j1];
}
// apply correction
sol[j1] += alpha;
sol[j2] -= alpha;
// update the local residual for corrections
for (size_t p = 0; p < nCorners; p++)
rhs[p] -= (A[p][j1] - A[p][j2]) * alpha;
}
}
}
// Point-to-Polyhedron
template <class Coordinate>
double distance(const Coordinate &target,
const ConvexPolyhedron<Coordinate> &segment,
double correctionTolerance) {
size_t nCorners = segment.vertices.size();
double const barycentricTolerance = 1e-14;
Dune::DynamicMatrix<double> A(nCorners, nCorners);
populateMatrix(segment, A);
Dune::DynamicVector<double> l(nCorners);
for (size_t i = 0; i < nCorners; i++)
l[i] = target * segment.vertices[i];
// choose a feasible initial solution
Dune::DynamicVector<double> sol(nCorners);
for (size_t i = 0; i < nCorners; i++)
sol[i] = 1.0 / nCorners;
Coordinate oldCoordinates = segment.barycentricToGlobal(sol);
Coordinate newCoordinates;
size_t maxIterations = 1000;
size_t iterations;
for (iterations = 0; iterations < maxIterations; ++iterations) {
approximate(target, A, l, segment, sol);
newCoordinates = segment.barycentricToGlobal(sol);
if ((oldCoordinates - newCoordinates).two_norm() < correctionTolerance)
break;
oldCoordinates = newCoordinates;
}
assert(iterations < maxIterations);
segment.sanityCheck(sol, barycentricTolerance);
return (target - newCoordinates).two_norm();
}
// Polyhedron-to-Polyhedron
template <class Coordinate>
double distance(const ConvexPolyhedron<Coordinate> &s1,
const ConvexPolyhedron<Coordinate> &s2,
double valueCorrectionTolerance) {
size_t nCorners1 = s1.vertices.size();
size_t nCorners2 = s2.vertices.size();
double const barycentricTolerance = 1e-14;
// choose feasible initial solutions
Dune::DynamicVector<double> sol1(nCorners1);
for (size_t i = 0; i < s1.vertices.size(); i++)
sol1[i] = 1.0 / s1.vertices.size();
auto c1 = s1.barycentricToGlobal(sol1);
Dune::DynamicVector<double> sol2(nCorners2);
for (size_t i = 0; i < nCorners2; i++)
sol2[i] = 1.0 / nCorners2;
auto c2 = s2.barycentricToGlobal(sol2);
double oldDistance = (c2 - c1).two_norm();
Dune::DynamicMatrix<double> A1(nCorners1, nCorners1);
populateMatrix(s1, A1);
Dune::DynamicMatrix<double> A2(nCorners2, nCorners2);
populateMatrix(s2, A2);
size_t maxIterations = 1000;
size_t iterations;
for (iterations = 0; iterations < maxIterations; ++iterations) {
Dune::DynamicVector<double> l2(nCorners2);
for (size_t i = 0; i < nCorners2; i++)
l2[i] = c1 * s2.vertices[i];
approximate(c1, A2, l2, s2, sol2);
c2 = s2.barycentricToGlobal(sol2);
Dune::DynamicVector<double> l1(nCorners1);
for (size_t i = 0; i < nCorners1; i++)
l1[i] = c2 * s1.vertices[i];
approximate(c2, A1, l1, s1, sol1);
c1 = s1.barycentricToGlobal(sol1);
double const newDistance = (c2 - c1).two_norm();
assert(newDistance - oldDistance <= 1e-12); // debugging
if (oldDistance - newDistance <= valueCorrectionTolerance) {
s1.sanityCheck(sol1, barycentricTolerance);
s2.sanityCheck(sol2, barycentricTolerance);
return newDistance;
}
oldDistance = newDistance;
}
assert(false);
}
#endif
dune/tectonic/quadraticenergy.hh
View file @ 8fe950bc
...@@ -8,12 +8,11 @@ template <class NonlinearityTEMPLATE> class QuadraticEnergy { ...@@ -8,12 +8,11 @@ template <class NonlinearityTEMPLATE> class QuadraticEnergy {
using Nonlinearity = NonlinearityTEMPLATE; using Nonlinearity = NonlinearityTEMPLATE;
using LocalVector = typename Nonlinearity::VectorType; using LocalVector = typename Nonlinearity::VectorType;
QuadraticEnergy(double alpha, LocalVector const &b, QuadraticEnergy(double alpha, LocalVector const &b, Nonlinearity const &phi)
std::shared_ptr<Nonlinearity const> phi)
: alpha(alpha), b(b), phi(phi) {} : alpha(alpha), b(b), phi(phi) {}
double const alpha; double const alpha;
LocalVector const &b; LocalVector const &b;
std::shared_ptr<Nonlinearity const> const phi; Nonlinearity const &phi;
}; };
#endif #endif
dune/tectonic/test/CMakeLists.txt deleted 100644 → 0
View file @ 4bea0810
set(TESTS test-polyhedral-minimisation)
foreach(_test ${TESTS})
add_executable(${_test} EXCLUDE_FROM_ALL ${_test}.cc)
target_link_dune_default_libraries(${_test})
add_test(${_test} ${_test})
endforeach()
add_directory_test_target(_test_target)
add_dependencies(${_test_target} ${TESTS})
dune/tectonic/test/test-polyhedral-minimisation.cc deleted 100644 → 0
View file @ 4bea0810
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/fvector.hh>
#include <dune/tectonic/polyhedrondistance.hh>
int main() {
using LocalVector = Dune::FieldVector<double, 2>;
auto const test =
[&](ConvexPolyhedron<LocalVector> const &p1,
ConvexPolyhedron<LocalVector> const &p2, double analyticalDistance) {
LocalVector target;
{
double const error =
std::abs(analyticalDistance - distance(p1, p2, 1e-12));
std::cout << "error: " << error << std::endl;
assert(error < 1e-12);
}
};
{
/*
* Calculate the distance between two triangles, where it is attained at a
* face-vertex pair
*
* O
* |\
* | \
* O O--O
* |\
* | \
* O--O
*/
double const analyticalDistance = std::sqrt(2.0);
LocalVector tmp;
ConvexPolyhedron<LocalVector> bs1;
bs1.vertices.resize(3);
tmp[0] = 0;
tmp[1] = 0;
bs1.vertices[0] = tmp;
tmp[0] = 0;
tmp[1] = 2;
bs1.vertices[1] = tmp;
tmp[0] = 2;
tmp[1] = 0;
bs1.vertices[2] = tmp;
ConvexPolyhedron<LocalVector> bs2;
bs2.vertices.resize(3);
tmp[0] = 2;
tmp[1] = 2;
bs2.vertices[0] = tmp;
tmp[0] = 2;
tmp[1] = 4;
bs2.vertices[1] = tmp;
tmp[0] = 4;
tmp[1] = 2;
bs2.vertices[2] = tmp;
test(bs1, bs2, analyticalDistance);
}
{
/*
* Calculate the distance between two triangles, where it is
* attained in a face-face pair
*
* O--O
* \ |
* \|
* O O
* |\
* | \
* O--O
*/
double const analyticalDistance = 2.0 * std::sqrt(2.0);
LocalVector tmp;
ConvexPolyhedron<LocalVector> bs1;
bs1.vertices.resize(3);
tmp[0] = 0;
tmp[1] = 0;
bs1.vertices[0] = tmp;
tmp[0] = 0;
tmp[1] = 2;
bs1.vertices[1] = tmp;
tmp[0] = 2;
tmp[1] = 0;
bs1.vertices[2] = tmp;
ConvexPolyhedron<LocalVector> bs2;
bs2.vertices.resize(3);
tmp[0] = 4;
tmp[1] = 4;
bs2.vertices[0] = tmp;
tmp[0] = 2;
tmp[1] = 4;
bs2.vertices[1] = tmp;
tmp[0] = 4;
tmp[1] = 2;
bs2.vertices[2] = tmp;
test(bs1, bs2, analyticalDistance);
}
{
/*
* Calculate the distance between two triangles, where it is
* attained in a vertex-vertex pair
*
* O
* |\
* | \
* O--O O--O
* \ |
* \|
* O
*/
double analyticalDistance = 2.0;
LocalVector tmp;
ConvexPolyhedron<LocalVector> bs1;
bs1.vertices.resize(3);
tmp[0] = 2;
tmp[1] = 2;
bs1.vertices[0] = tmp;
tmp[0] = 0;
tmp[1] = 2;
bs1.vertices[1] = tmp;
tmp[0] = 2;
tmp[1] = 0;
bs1.vertices[2] = tmp;
ConvexPolyhedron<LocalVector> bs2;
bs2.vertices.resize(3);
tmp[0] = 4;
tmp[1] = 2;
bs2.vertices[0] = tmp;
tmp[0] = 4;
tmp[1] = 4;
bs2.vertices[1] = tmp;
tmp[0] = 6;
tmp[1] = 2;
bs2.vertices[2] = tmp;
test(bs1, bs2, analyticalDistance);
}
}
src/CMakeLists.txt
View file @ 8fe950bc
set(SOURCE_FILES set(SW_SOURCE_FILES
assemblers.cc assemblers.cc
boundary_writer.cc
coupledtimestepper.cc
enumparser.cc enumparser.cc
fixedpointiterator.cc hdf5/frictionalboundary-writer.cc
friction_writer.cc hdf5/iteration-writer.cc
sand-wedge.cc hdf5/patchinfo-writer.cc
sand-wedge-data/mygeometry.cc hdf5/restart-io.cc
sand-wedge-data/mygrid.cc hdf5/surface-writer.cc
solverfactory.cc hdf5/time-writer.cc
state.cc one-body-problem-data/mygeometry.cc
timestepping.cc one-body-problem-data/mygrid.cc
one-body-problem.cc
spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc
time-stepping/adaptivetimestepper.cc
time-stepping/coupledtimestepper.cc
time-stepping/rate.cc
time-stepping/rate/rateupdater.cc
time-stepping/state.cc
vtk.cc
)
set(UGW_SOURCE_FILES
assemblers.cc # FIXME
one-body-problem-data/mygrid.cc
uniform-grid-writer.cc
vtk.cc vtk.cc
) )
file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/sand-wedge-data")
dune_symlink_to_source_files("sand-wedge-data/boundaryconditions.py")
dune_symlink_to_source_files("sand-wedge-data/parset.cfg")
find_package(Boost REQUIRED system filesystem)
include_directories(${Boost_INCLUDE_DIR})
foreach(_dim 2 3) foreach(_dim 2 3)
set(_target sand-wedge-${_dim}D) set(_sw_target one-body-problem-${_dim}D)
add_executable(${_target} ${SOURCE_FILES}) set(_ugw_target uniform-grid-writer-${_dim}D)
add_dune_pythonlibs_flags(${_target})
add_dune_ug_flags(${_target}) add_executable(${_sw_target} ${SW_SOURCE_FILES})
add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
add_dune_ug_flags(${_sw_target})
add_dune_ug_flags(${_ugw_target})
target_link_libraries(${_target} ${Boost_FILESYSTEM_LIBRARY}) add_dune_hdf5_flags(${_sw_target})
target_link_libraries(${_target} ${Boost_SYSTEM_LIBRARY})
set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}") set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
endforeach() endforeach()
src/adaptivetimestepper.hh deleted 100644 → 0
View file @ 4bea0810
#include "coupledtimestepper.hh"
template <typename T1, typename T2>
std::pair<T1, T2> clonePair(std::pair<T1, T2> in) {
return { in.first->clone(), in.second->clone() };
}
template <class Factory, class UpdaterPair, class ErrorNorm>
class AdaptiveTimeStepper {
using StateUpdater = typename UpdaterPair::first_type::element_type;
using VelocityUpdater = typename UpdaterPair::second_type::element_type;
using Vector = typename Factory::Vector;
using ConvexProblem = typename Factory::ConvexProblem;
using Nonlinearity = typename ConvexProblem::NonlinearityType;
using MyCoupledTimeStepper =
CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>;
public:
AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, UpdaterPair &current,
std::function<void(double, Vector &)> externalForces,
ErrorNorm const &errorNorm,
std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine)
: finalTime_(parset.get<double>("problem.finalTime")),
relativeTime_(0.0),
relativeTau_(1e-6), // FIXME (not really important, though)
factory_(factory),
parset_(parset),
globalFriction_(globalFriction),
current_(current),
R1_(clonePair(current_)),
externalForces_(externalForces),
mustRefine_(mustRefine),
iterationWriter_("iterations", std::fstream::out),
errorNorm_(errorNorm) {
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, R1_.first, R1_.second,
errorNorm, externalForces_);
stepAndReport("R1", coupledTimeStepper, relativeTime_, relativeTau_);
iterationWriter_ << std::endl;
}
bool reachedEnd() { return relativeTime_ >= 1.0; }
bool coarsen() {
bool didCoarsen = false;
while (relativeTime_ + relativeTau_ < 1.0) {
R2_ = clonePair(R1_);
{
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, R2_.first,
R2_.second, errorNorm_, externalForces_);
stepAndReport("R2", coupledTimeStepper, relativeTime_ + relativeTau_,
relativeTau_);
}
UpdaterPair C = clonePair(current_);
{
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, C.first, C.second,
errorNorm_, externalForces_);
stepAndReport("C", coupledTimeStepper, relativeTime_,
2.0 * relativeTau_);
}
if (!mustRefine_(C, R2_)) {
R2_ = { nullptr, nullptr };
R1_ = C;
relativeTau_ *= 2.0;
didCoarsen = true;
} else {
break;
}
}
return didCoarsen;
}
void refine() {
while (true) {
UpdaterPair F2 = clonePair(current_);
UpdaterPair F1;
{
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, F2.first, F2.second,
errorNorm_, externalForces_);
stepAndReport("F1", coupledTimeStepper, relativeTime_,
relativeTau_ / 2.0);
F1 = clonePair(F2);
stepAndReport("F2", coupledTimeStepper,
relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
}
if (!mustRefine_(R1_, F2)) {
break;
} else {
R1_ = F1;
R2_ = F2;
relativeTau_ /= 2.0;
}
}
}
void advance() {
if (!coarsen())
refine();
iterationWriter_ << std::endl;
current_ = R1_;
R1_ = R2_;
relativeTime_ += relativeTau_;
}
double getRelativeTime() { return relativeTime_; }
double getRelativeTau() { return relativeTau_; }
private:
void stepAndReport(std::string type, MyCoupledTimeStepper &stepper,
double rTime, double rTau) {
iterationWriter_ << type << " " << stepper.step(rTime, rTau) << " "
<< std::flush;
}
double finalTime_;
double relativeTime_;
double relativeTau_;
Factory &factory_;
Dune::ParameterTree const &parset_;
std::shared_ptr<Nonlinearity> globalFriction_;
UpdaterPair &current_;
UpdaterPair R1_;
UpdaterPair R2_;
std::function<void(double, Vector &)> externalForces_;
std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine_;
std::fstream iterationWriter_;
ErrorNorm const &errorNorm_;
};
src/assemblers.cc
View file @ 8fe950bc
...@@ -7,14 +7,15 @@ ...@@ -7,14 +7,15 @@
#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh> #include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh> #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh> #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh> #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh> #include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh> #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh> #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
#include <dune/fufem/boundarypatch.hh> #include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/computestress.hh>
#include <dune/fufem/functions/basisgridfunction.hh> #include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/constantfunction.hh> #include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functiontools/p0p1interpolation.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh> #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/tectonic/frictionpotential.hh> #include <dune/tectonic/frictionpotential.hh>
...@@ -33,7 +34,7 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView) ...@@ -33,7 +34,7 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
template <class GridView, int dimension> template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass( void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
BoundaryPatch<GridView> const &frictionalBoundary, BoundaryPatch<GridView> const &frictionalBoundary,
ScalarMatrix &frictionalBoundaryMass) { ScalarMatrix &frictionalBoundaryMass) const {
BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis, BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
frictionalBoundaryMassAssembler(frictionalBoundary); frictionalBoundaryMassAssembler(frictionalBoundary);
...@@ -45,7 +46,7 @@ template <class GridView, int dimension> ...@@ -45,7 +46,7 @@ template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleMass( void MyAssembler<GridView, dimension>::assembleMass(
Dune::VirtualFunction<LocalVector, LocalScalarVector> const & Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
densityFunction, densityFunction,
Matrix &M) { Matrix &M) const {
// NOTE: We treat the weight as a constant function // NOTE: We treat the weight as a constant function
QuadratureRuleKey quadKey(dimension, 0); QuadratureRuleKey quadKey(dimension, 0);
...@@ -58,7 +59,7 @@ void MyAssembler<GridView, dimension>::assembleMass( ...@@ -58,7 +59,7 @@ void MyAssembler<GridView, dimension>::assembleMass(
template <class GridView, int dimension> template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu, void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
Matrix &A) { Matrix &A) const {
StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
localStiffness(E, nu); localStiffness(E, nu);
vertexAssembler.assembleOperator(localStiffness, A); vertexAssembler.assembleOperator(localStiffness, A);
...@@ -68,7 +69,7 @@ template <class GridView, int dimension> ...@@ -68,7 +69,7 @@ template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleViscosity( void MyAssembler<GridView, dimension>::assembleViscosity(
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity, Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity, Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
Matrix &C) { Matrix &C) const {
// NOTE: We treat the weights as constant functions // NOTE: We treat the weights as constant functions
QuadratureRuleKey shearViscosityKey(dimension, 0); QuadratureRuleKey shearViscosityKey(dimension, 0);
QuadratureRuleKey bulkViscosityKey(dimension, 0); QuadratureRuleKey bulkViscosityKey(dimension, 0);
...@@ -83,7 +84,7 @@ void MyAssembler<GridView, dimension>::assembleViscosity( ...@@ -83,7 +84,7 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
template <class GridView, int dimension> template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleBodyForce( void MyAssembler<GridView, dimension>::assembleBodyForce(
Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField, Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
Vector &f) { Vector &f) const {
L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector> L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
gravityFunctionalAssembler(gravityField); gravityFunctionalAssembler(gravityField);
vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f); vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
...@@ -92,7 +93,8 @@ void MyAssembler<GridView, dimension>::assembleBodyForce( ...@@ -92,7 +93,8 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
template <class GridView, int dimension> template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNeumann( void MyAssembler<GridView, dimension>::assembleNeumann(
BoundaryPatch<GridView> const &neumannBoundary, Vector &f, BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
Dune::VirtualFunction<double, double> const &neumann, double relativeTime) { Dune::VirtualFunction<double, double> const &neumann,
double relativeTime) const {
LocalVector localNeumann(0); LocalVector localNeumann(0);
neumann.evaluate(relativeTime, localNeumann[0]); neumann.evaluate(relativeTime, localNeumann[0]);
NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler( NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
...@@ -103,25 +105,35 @@ void MyAssembler<GridView, dimension>::assembleNeumann( ...@@ -103,25 +105,35 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
} }
template <class GridView, int dimension> template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNormalStress( void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
BoundaryPatch<GridView> const &frictionalBoundary, BoundaryPatch<GridView> const &frictionalBoundary,
ScalarVector &normalStress, double youngModulus, double poissonRatio, ScalarVector &weightedNormalStress, double youngModulus,
Vector const &displacement) { double poissonRatio, Vector const &displacement) const {
Vector traction; BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!) displacement);
(frictionalBoundary, displacement, traction, youngModulus, poissonRatio); Vector traction(cellBasis.size());
NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
std::vector<typename Vector::block_type> normals; youngModulus, poissonRatio, &displacementFunction, 1);
frictionalBoundary.getNormals(normals); cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
for (size_t i = 0; i < traction.size(); ++i) { frictionalBoundary);
normalStress[i] = normals[i] * traction[i];
if (normalStress[i] > 0.0) { auto const nodalTractionAverage =
normalStress[i] = 0.0; interpolateP0ToP1(frictionalBoundary, traction);
std::cout << "Warning: Manually reducing positive normal stress to zero."
<< std::endl; ScalarVector weights;
} {
NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
frictionalBoundaryAssembler(
std::make_shared<ConstantFunction<
LocalVector, typename ScalarVector::block_type>>(1));
vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
weights, frictionalBoundary);
} }
auto const normals = frictionalBoundary.getNormals();
for (size_t i = 0; i < vertexBasis.size(); ++i)
weightedNormalStress[i] =
std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
} }
template <class GridView, int dimension> template <class GridView, int dimension>
...@@ -129,9 +141,9 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( ...@@ -129,9 +141,9 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
Config::FrictionModel frictionModel, Config::FrictionModel frictionModel,
BoundaryPatch<GridView> const &frictionalBoundary, BoundaryPatch<GridView> const &frictionalBoundary,
GlobalFrictionData<dimension> const &frictionInfo, GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &normalStress) ScalarVector const &weightedNormalStress) const
-> std::shared_ptr<GlobalFriction<Matrix, Vector>> { -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
// Lump negative normal stress (kludge) // Lumping of the nonlinearity
ScalarVector weights; ScalarVector weights;
{ {
NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type> NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
...@@ -145,11 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( ...@@ -145,11 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
case Config::Truncated: case Config::Truncated:
return std::make_shared<GlobalRateStateFriction< return std::make_shared<GlobalRateStateFriction<
Matrix, Vector, TruncatedRateState, GridView>>( Matrix, Vector, TruncatedRateState, GridView>>(
frictionalBoundary, gridView, frictionInfo, weights, normalStress); frictionalBoundary, frictionInfo, weights, weightedNormalStress);
case Config::Regularised: case Config::Regularised:
return std::make_shared<GlobalRateStateFriction< return std::make_shared<GlobalRateStateFriction<
Matrix, Vector, RegularisedRateState, GridView>>( Matrix, Vector, RegularisedRateState, GridView>>(
frictionalBoundary, gridView, frictionInfo, weights, normalStress); frictionalBoundary, frictionInfo, weights, weightedNormalStress);
default: default:
assert(false); assert(false);
} }
...@@ -158,7 +170,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( ...@@ -158,7 +170,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
template <class GridView, int dimension> template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleVonMisesStress( void MyAssembler<GridView, dimension>::assembleVonMisesStress(
double youngModulus, double poissonRatio, Vector const &u, double youngModulus, double poissonRatio, Vector const &u,
ScalarVector &stress) { ScalarVector &stress) const {
auto const gridDisplacement = auto const gridDisplacement =
std::make_shared<BasisGridFunction<VertexBasis, Vector> const>( std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
vertexBasis, u); vertexBasis, u);
......