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 (63)
Hide whitespace changes
Inline Side-by-side
Showing
with 229 additions and 703 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} 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
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
...@@ -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);
} }
...@@ -102,18 +104,18 @@ class RegularisedRateState : public FrictionPotential { ...@@ -102,18 +104,18 @@ class RegularisedRateState : public FrictionPotential {
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) const { 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,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 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], 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>
...@@ -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/pointtractionboundaryassembler.hh deleted 100644 → 0
View file @ 549f8113
#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/polyhedrondistance.hh deleted 100644 → 0
View file @ 549f8113
#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 @ 549f8113
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 @ 549f8113
#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
adaptivetimestepper.cc
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
rate.cc hdf5/patchinfo-writer.cc
rate/rateupdater.cc hdf5/restart-io.cc
sand-wedge.cc hdf5/surface-writer.cc
sand-wedge-data/mygeometry.cc hdf5/time-writer.cc
sand-wedge-data/mygrid.cc one-body-problem-data/mygeometry.cc
solverfactory.cc one-body-problem-data/mygrid.cc
state.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 serialization)
include_directories(${Boost_INCLUDE_DIR})
# dataio.hh expects this to be set
set_property(DIRECTORY APPEND PROPERTY COMPILE_DEFINITIONS "HAVE_BOOST_SERIALIZATION")
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_SERIALIZATION_LIBRARY})
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/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> weightedNormalStressBoundaryAssembler(
youngModulus, poissonRatio, &displacementFunction, 1); youngModulus, poissonRatio, &displacementFunction, 1);
vertexAssembler.assembleBoundaryFunctional( cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
weightedNormalStressBoundaryAssembler, 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);
} }
......