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 (37)
Hide whitespace changes
Inline Side-by-side
Showing
with 189 additions and 297 deletions
CMakeLists.txt
View file @ 8fe950bc
...@@ -18,8 +18,7 @@ include(DuneMacros) ...@@ -18,8 +18,7 @@ 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(Boost REQUIRED system filesystem) find_package(HDF5 COMPONENTS C REQUIRED)
find_package(HDF5 REQUIRED)
add_subdirectory("src") add_subdirectory("src")
add_subdirectory("dune") add_subdirectory("dune")
......
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} Libs: -L${libdir}
Cflags: -I${includedir} Cflags: -I${includedir}
dune.module
View file @ 8fe950bc
...@@ -3,6 +3,6 @@ ...@@ -3,6 +3,6 @@
################################ ################################
Module: dune-tectonic Module: dune-tectonic
Version: 2.4-1 Version: 2.5-1
Maintainer: elias.pipping@fu-berlin.de Maintainer: elias.pipping@fu-berlin.de
Depends: dune-common dune-fufem dune-geometry dune-grid dune-istl dune-solvers dune-tnnmg Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
dune/tectonic/CMakeLists.txt
View file @ 8fe950bc
...@@ -10,7 +10,6 @@ install(FILES ...@@ -10,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
...@@ -33,32 +33,32 @@ class TruncatedRateState : public FrictionPotential { ...@@ -33,32 +33,32 @@ class TruncatedRateState : public FrictionPotential {
FrictionData _fd) FrictionData _fd)
: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {} : 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 - weightedNormalStress * 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 -weightedNormalStress * (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);
} }
...@@ -76,21 +76,23 @@ class RegularisedRateState : public FrictionPotential { ...@@ -76,21 +76,23 @@ class RegularisedRateState : public FrictionPotential {
FrictionData _fd) FrictionData _fd)
: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {} : 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 - weightedNormalStress * coefficientOfFriction(V); return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
} }
double second_deriv(double V) const { double second_deriv(double V) const override {
return -weightedNormalStress * 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);
} }
...@@ -104,16 +106,16 @@ class RegularisedRateState : public FrictionPotential { ...@@ -104,16 +106,16 @@ class RegularisedRateState : public FrictionPotential {
class ZeroFunction : 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,18 @@ template <class Matrix, class Vector> class GlobalFriction { ...@@ -71,21 +66,18 @@ 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);
} }
ScalarVector coefficientOfFriction(Vector const &x) const { ScalarVector coefficientOfFriction(Vector const &x) const {
ScalarVector ret(x.size()); ScalarVector ret(x.size());
for (size_t i = 0; i < x.size(); ++i) for (size_t i = 0; i < x.size(); ++i)
ret[i] = restriction(i)->coefficientOfFriction(x[i]); ret[i] = restriction(i).coefficientOfFriction(x[i]);
return ret; return ret;
} }
......
dune/tectonic/globalratestatefriction.hh
View file @ 8fe950bc
...@@ -10,57 +10,62 @@ ...@@ -10,57 +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,
ScalarVector const &weights, ScalarVector const &weights,
ScalarVector const &weightedNormalStress) ScalarVector const &weightedNormalStress)
: restrictions(weightedNormalStress.size()) { : restrictions(), localToGlobal(), zeroFriction() {
auto zeroNonlinearity = auto const gridView = frictionalBoundary.gridView();
std::make_shared<Friction>(std::make_shared<ZeroFunction>());
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] = zeroNonlinearity;
continue; continue;
}
auto const fp = std::make_shared<ScalarFriction>( localToGlobal.emplace_back(i);
weights[i], weightedNormalStress[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>
...@@ -256,7 +255,8 @@ class MyBlockProblem<ConvexProblem>::IterateObject { ...@@ -256,7 +255,8 @@ 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_);
} }
......
dune/tectonic/pointtractionboundaryassembler.hh deleted 100644 → 0
View file @ 8c878ee4
#ifndef DUNE_TECTONIC_POINTTRACTIONBOUNDARYASSEMBLER_HH
#define DUNE_TECTONIC_POINTTRACTIONBOUNDARYASSEMBLER_HH
// Based on
// dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh
#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/assemblers/localboundaryassembler.hh>
#include <dune/fufem/symmetrictensor.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
/** \brief Computes the surface traction of a linear elastic material,
pointwise. This is done in such a manner that the sum over all
vertices equals the (approximate) integral */
template <class GridType>
class PointTractionBoundaryAssembler
: public LocalBoundaryAssembler<
GridType,
Dune::FieldVector<typename GridType::ctype, GridType::dimension>> {
static const int dim = GridType::dimension;
typedef typename GridType::ctype ctype;
typedef typename Dune::FieldVector<ctype, dim> T;
public:
typedef typename LocalBoundaryAssembler<GridType, T>::LocalVector LocalVector;
typedef VirtualGridFunction<
GridType, Dune::FieldVector<ctype, GridType::dimension>> GridFunction;
//! Create assembler for grid
PointTractionBoundaryAssembler(double E, double nu,
const GridFunction* displacement, int order)
: E_(E), nu_(nu), displacement_(displacement), order_(order) {}
/** \brief Assemble the normal stress at a boundary face.*/
template <class TrialLocalFE, class BoundaryIterator>
void assemble(const BoundaryIterator& it, LocalVector& localVector,
const TrialLocalFE& tFE) {
localVector = 0.0;
// geometry of the boundary face
const typename BoundaryIterator::Intersection::Geometry segmentGeometry =
it->geometry();
// get quadrature rule
const int order = (it->inside()->type().isSimplex())
? 2 * (order_ - 1)
: 2 * (order_ * dim - 1);
// get quadrature rule
const Dune::QuadratureRule<ctype, dim - 1>& quad =
QuadratureRuleCache<ctype, dim - 1>::rule(segmentGeometry.type(), order,
false);
// loop over quadrature points
for (size_t pt = 0; pt < quad.size(); ++pt) {
// get quadrature point
const Dune::FieldVector<ctype, dim - 1>& quadPos = quad[pt].position();
// get integration factor
const ctype integrationElement =
segmentGeometry.integrationElement(quadPos);
// get outer unit normal at quad pos
Dune::FieldVector<ctype, dim> normalAtQuadPos =
it->unitOuterNormal(quadPos);
// position of the quadrature point within the element
const Dune::FieldVector<ctype, dim> elementQuadPos =
it->geometryInInside().global(quadPos);
// evaluate the displacement gradient at quad point of the element
typename GridFunction::DerivativeType localDispGradient;
if (displacement_->isDefinedOn(*it->inside()))
displacement_->evaluateDerivativeLocal(*it->inside(), elementQuadPos,
localDispGradient);
else
displacement_->evaluateDerivative(segmentGeometry.global(quadPos),
localDispGradient);
// Compute strain tensor from the displacement gradient
SymmetricTensor<dim> strain;
computeStrain(localDispGradient, strain);
// and the stress
SymmetricTensor<dim> stress = hookeTimesStrain(strain);
normalAtQuadPos *= quad[pt].weight() * integrationElement;
typename LocalVector::block_type tractionAtQuadPos;
stress.umv(normalAtQuadPos, tractionAtQuadPos);
// Divide the traction
tractionAtQuadPos /= segmentGeometry.corners();
for (size_t k = 0; k < localVector.size(); k++)
localVector[k] += tractionAtQuadPos;
}
}
~PointTractionBoundaryAssembler() {}
private:
/** \brief Young's modulus */
double E_;
/** \brief The Poisson ratio */
double nu_;
/** \brief The displacement.*/
const GridFunction* displacement_;
/** \brief The gridfunction is usually corresponding to some finite element
* solution. This is its order.*/
int order_;
/** \brief Compute linear strain tensor from the deformation gradient. */
void computeStrain(const Dune::FieldMatrix<double, dim, dim>& grad,
SymmetricTensor<dim>& strain) const {
for (int i = 0; i < dim; ++i) {
strain(i, i) = grad[i][i];
for (int j = i + 1; j < dim; ++j)
strain(i, j) = 0.5 * (grad[i][j] + grad[j][i]);
}
}
/** \brief Compute the stress tensor from strain tensor. */
SymmetricTensor<dim> hookeTimesStrain(
const SymmetricTensor<dim>& strain) const {
SymmetricTensor<dim> stress = strain;
stress *= E_ / (1 + nu_);
double f = E_ * nu_ / ((1 + nu_) * (1 - 2 * nu_)) * strain.trace();
for (int i = 0; i < dim; i++)
stress[i] += f;
return stress;
}
};
#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
src/CMakeLists.txt
View file @ 8fe950bc
set(SOURCE_FILES set(SW_SOURCE_FILES
assemblers.cc assemblers.cc
enumparser.cc enumparser.cc
hdf5/frictionalboundary-writer.cc hdf5/frictionalboundary-writer.cc
...@@ -7,9 +7,9 @@ set(SOURCE_FILES ...@@ -7,9 +7,9 @@ set(SOURCE_FILES
hdf5/restart-io.cc hdf5/restart-io.cc
hdf5/surface-writer.cc hdf5/surface-writer.cc
hdf5/time-writer.cc hdf5/time-writer.cc
sand-wedge-data/mygeometry.cc one-body-problem-data/mygeometry.cc
sand-wedge-data/mygrid.cc one-body-problem-data/mygrid.cc
sand-wedge.cc one-body-problem.cc
spatial-solving/fixedpointiterator.cc spatial-solving/fixedpointiterator.cc
spatial-solving/solverfactory.cc spatial-solving/solverfactory.cc
time-stepping/adaptivetimestepper.cc time-stepping/adaptivetimestepper.cc
...@@ -19,30 +19,24 @@ set(SOURCE_FILES ...@@ -19,30 +19,24 @@ set(SOURCE_FILES
time-stepping/state.cc time-stepping/state.cc
vtk.cc vtk.cc
) )
set(UGW_SOURCE_FILES
file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/sand-wedge-data") assemblers.cc # FIXME
dune_symlink_to_source_files("sand-wedge-data/boundaryconditions.py") one-body-problem-data/mygrid.cc
dune_symlink_to_source_files("sand-wedge-data/parset.cfg") uniform-grid-writer.cc
dune_symlink_to_source_files("sand-wedge-data/parset-2D.cfg") vtk.cc
dune_symlink_to_source_files("sand-wedge-data/parset-3D.cfg") )
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})
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS}) 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 INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIR}) set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
target_link_libraries(${_target} ${HDF5_LIBRARIES}) set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
target_link_libraries(${_target} "dl") # FIXME: missing from HDF5_LIBRARIES
set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
endforeach() endforeach()
# Mark UG as required
find_package(UG REQUIRED)
src/assemblers.cc
View file @ 8fe950bc
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
#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>
...@@ -14,11 +15,11 @@ ...@@ -14,11 +15,11 @@
#include <dune/fufem/boundarypatch.hh> #include <dune/fufem/boundarypatch.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>
#include <dune/tectonic/globalratestatefriction.hh> #include <dune/tectonic/globalratestatefriction.hh>
#include <dune/tectonic/pointtractionboundaryassembler.hh>
#include "assemblers.hh" #include "assemblers.hh"
...@@ -111,24 +112,28 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress( ...@@ -111,24 +112,28 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis, BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
displacement); displacement);
Vector traction; Vector traction(cellBasis.size());
// Note: We assume v = 0 here so that there is no viscous NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
// contribution to the stress.
PointTractionBoundaryAssembler<Grid> weightedTractionBoundaryAssembler(
youngModulus, poissonRatio, &displacementFunction, 1); youngModulus, poissonRatio, &displacementFunction, 1);
vertexAssembler.assembleBoundaryFunctional(weightedTractionBoundaryAssembler, cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
traction, frictionalBoundary); frictionalBoundary);
std::vector<typename Vector::block_type> normals; auto const nodalTractionAverage =
frictionalBoundary.getNormals(normals); interpolateP0ToP1(frictionalBoundary, traction);
for (size_t i = 0; i < traction.size(); ++i) {
weightedNormalStress[i] = normals[i] * traction[i]; ScalarVector weights;
if (weightedNormalStress[i] > 0.0) { {
weightedNormalStress[i] = 0.0; NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
std::cout << "Warning: Manually reducing positive normal stress to zero." frictionalBoundaryAssembler(
<< std::endl; 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>
...@@ -152,13 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity( ...@@ -152,13 +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, frictionalBoundary, frictionInfo, weights, weightedNormalStress);
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, frictionalBoundary, frictionInfo, weights, weightedNormalStress);
weightedNormalStress);
default: default:
assert(false); assert(false);
} }
......
src/enumparser.cc
View file @ 8fe950bc
...@@ -9,8 +9,8 @@ ...@@ -9,8 +9,8 @@
#include "enumparser.hh" #include "enumparser.hh"
template <class Enum> template <class Enum>
typename Dune::enable_if< typename std::enable_if<
!Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value, !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
std::istream &>::type std::istream &>::type
operator>>(std::istream &lhs, Enum &e) { operator>>(std::istream &lhs, Enum &e) {
std::string s; std::string s;
......
src/enumparser.hh
View file @ 8fe950bc
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// Copyright Carsten Graeser 2012 // Copyright Carsten Graeser 2012
#include <dune/common/typetraits.hh> #include <type_traits>
#include <dune/solvers/solvers/solver.hh> #include <dune/solvers/solvers/solver.hh>
...@@ -28,8 +28,8 @@ template <> struct StringToEnum<Config::PatchType> { ...@@ -28,8 +28,8 @@ template <> struct StringToEnum<Config::PatchType> {
}; };
template <class Enum> template <class Enum>
typename Dune::enable_if< typename std::enable_if<
!Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value, !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
std::istream &>::type std::istream &>::type
operator>>(std::istream &lhs, Enum &e); operator>>(std::istream &lhs, Enum &e);
#endif #endif
src/gridselector.hh
View file @ 8fe950bc
...@@ -12,10 +12,7 @@ ...@@ -12,10 +12,7 @@
#if !HAVE_ALUGRID #if !HAVE_ALUGRID
#error ALUGRID was requested but not found #error ALUGRID was requested but not found
#endif #endif
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wignored-qualifiers"
#include <dune/grid/alugrid.hh> #include <dune/grid/alugrid.hh>
#pragma clang diagnostic pop
using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>; using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
#elif WANT_GRID == WANT_UG #elif WANT_GRID == WANT_UG
......
src/hdf5-writer.hh
View file @ 8fe950bc
...@@ -3,12 +3,11 @@ ...@@ -3,12 +3,11 @@
#include <dune/fufem/functions/basisgridfunction.hh> #include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh> #include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/hdf5/hdf5file.hh> #include <dune/fufem/hdf5/file.hh>
#include "hdf5/frictionalboundary-writer.hh" #include "hdf5/frictionalboundary-writer.hh"
#include "hdf5/iteration-writer.hh" #include "hdf5/iteration-writer.hh"
#include "hdf5/patchinfo-writer.hh" #include "hdf5/patchinfo-writer.hh"
#include "hdf5/restart-io.hh"
#include "hdf5/surface-writer.hh" #include "hdf5/surface-writer.hh"
#include "hdf5/time-writer.hh" #include "hdf5/time-writer.hh"
...@@ -21,7 +20,7 @@ class HDF5Writer { ...@@ -21,7 +20,7 @@ class HDF5Writer {
using LocalVector = typename Vector::block_type; using LocalVector = typename Vector::block_type;
public: public:
HDF5Writer(HDF5Grouplike &file, Vector const &vertexCoordinates, HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates,
VertexBasis const &vertexBasis, Patch const &surface, VertexBasis const &vertexBasis, Patch const &surface,
Patch const &frictionalBoundary, Patch const &frictionalBoundary,
ConvexPolyhedron<LocalVector> const &weakPatch) ConvexPolyhedron<LocalVector> const &weakPatch)
...@@ -32,8 +31,8 @@ class HDF5Writer { ...@@ -32,8 +31,8 @@ class HDF5Writer {
patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch), patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
#endif #endif
surfaceWriter_(file_, vertexCoordinates, surface), surfaceWriter_(file_, vertexCoordinates, surface),
frictionalBoundaryWriter_(file_, vertexCoordinates, frictionalBoundary), frictionalBoundaryWriter_(file_, vertexCoordinates,
restartWriter_(file_, vertexCoordinates.size()) { frictionalBoundary) {
} }
template <class Friction> template <class Friction>
...@@ -53,12 +52,8 @@ class HDF5Writer { ...@@ -53,12 +52,8 @@ class HDF5Writer {
iterationWriter_.write(programState.timeStep, iterationCount); iterationWriter_.write(programState.timeStep, iterationCount);
} }
void writeRestart(ProgramState const &programState) {
restartWriter_.write(programState);
}
private: private:
HDF5Grouplike &file_; HDF5::Grouplike &file_;
IterationWriter iterationWriter_; IterationWriter iterationWriter_;
TimeWriter<ProgramState> timeWriter_; TimeWriter<ProgramState> timeWriter_;
...@@ -67,7 +62,5 @@ class HDF5Writer { ...@@ -67,7 +62,5 @@ class HDF5Writer {
#endif #endif
SurfaceWriter<ProgramState, GridView> surfaceWriter_; SurfaceWriter<ProgramState, GridView> surfaceWriter_;
FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_; FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
RestartIO<ProgramState> restartWriter_;
}; };
#endif #endif