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)
# start a dune project with information from dune.module
dune_project()
find_package(HDF5 COMPONENTS C REQUIRED)
add_subdirectory("src")
add_subdirectory("dune")
add_subdirectory("doc")
......
dune-tectonic.pc.in
View file @ 8fe950bc
......@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
Version: @VERSION@
Description: dune-tectonic module
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}
Cflags: -I${includedir}
dune.module
View file @ 8fe950bc
......@@ -3,6 +3,6 @@
################################
Module: dune-tectonic
Version: 2.4-1
Version: 2.5-1
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
body.hh
frictiondata.hh
......@@ -12,7 +10,6 @@ install(FILES
minimisation.hh
myblockproblem.hh
mydirectionalconvexfunction.hh
pointtractionboundaryassembler.hh
quadraticenergy.hh
tectonic.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
dune/tectonic/frictionpotential.hh
View file @ 8fe950bc
......@@ -33,32 +33,32 @@ class TruncatedRateState : public FrictionPotential {
FrictionData _fd)
: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
double coefficientOfFriction(double V) const {
double coefficientOfFriction(double V) const override {
if (V <= Vmin)
return 0.0;
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);
}
double second_deriv(double V) const {
double second_deriv(double V) const override {
if (V <= Vmin)
return 0;
return -weightedNormalStress * (fd.a / V);
}
double regularity(double V) const {
double regularity(double V) const override {
if (std::abs(V - Vmin) < 1e-14) // TODO
return std::numeric_limits<double>::infinity();
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;
Vmin = fd.V0 / std::exp(logrest);
}
......@@ -76,21 +76,23 @@ class RegularisedRateState : public FrictionPotential {
FrictionData _fd)
: 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);
}
double differential(double V) const {
double differential(double V) const override {
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);
}
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;
Vmin = fd.V0 / std::exp(logrest);
}
......@@ -102,18 +104,18 @@ class RegularisedRateState : public FrictionPotential {
double Vmin;
};
class TrivialFunction : public FrictionPotential {
class ZeroFunction : public FrictionPotential {
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
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 {
using LocalMatrix = typename Matrix::block_type;
using LocalVectorType = typename Vector::block_type;
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 tmp = 0;
for (size_t i = 0; i < x.size(); ++i) {
auto const res = restriction(i);
tmp += (*res)(x[i]);
tmp += restriction(i)(x[i]);
}
return tmp;
}
......@@ -36,14 +35,11 @@ template <class Matrix, class Vector> class GlobalFriction {
/*
Return a restriction of the outer function to the i'th node.
*/
std::shared_ptr<LocalFriction<block_size>> virtual restriction(
size_t i) const = 0;
LocalNonlinearity const virtual &restriction(size_t i) const = 0;
void addHessian(Vector const &v, Matrix &hessian) const {
for (size_t i = 0; i < v.size(); ++i) {
auto const res = restriction(i);
res->addHessian(v[i], hessian[i][i]);
}
for (size_t i = 0; i < v.size(); ++i)
restriction(i).addHessian(v[i], hessian[i][i]);
}
void directionalDomain(Vector const &, Vector const &,
......@@ -58,8 +54,7 @@ template <class Matrix, class Vector> class GlobalFriction {
subdifferential[0] = subdifferential[1] = 0;
for (size_t i = 0; i < u.size(); ++i) {
Dune::Solvers::Interval<double> D;
auto const res = restriction(i);
res->directionalSubDiff(u[i], v[i], D);
restriction(i).directionalSubDiff(u[i], v[i], D);
subdifferential[0] += D[0];
subdifferential[1] += D[1];
}
......@@ -71,21 +66,19 @@ template <class Matrix, class Vector> class GlobalFriction {
}
void addGradient(Vector const &v, Vector &gradient) const {
for (size_t i = 0; i < v.size(); ++i) {
auto const res = restriction(i);
res->addGradient(v[i], gradient[i]);
}
for (size_t i = 0; i < v.size(); ++i)
restriction(i).addGradient(v[i], gradient[i]);
}
double regularity(size_t i, typename Vector::block_type const &x) const {
auto const res = restriction(i);
return res->regularity(x);
return restriction(i).regularity(x);
}
void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) const {
coefficient.resize(x.size());
ScalarVector coefficientOfFriction(Vector const &x) const {
ScalarVector ret(x.size());
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;
......
dune/tectonic/globalratestatefriction.hh
View file @ 8fe950bc
......@@ -10,57 +10,62 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfrictiondata.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/index-in-sorted-range.hh>
template <class Matrix, class Vector, class ScalarFriction, class GridView>
class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
public:
using GlobalFriction<Matrix, Vector>::block_size;
using typename GlobalFriction<Matrix, Vector>::Friction;
using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
private:
using typename GlobalFriction<Matrix, Vector>::ScalarVector;
public:
GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
GridView const &gridView,
GlobalFrictionData<block_size> const &frictionInfo,
ScalarVector const &weights,
ScalarVector const &weightedNormalStress)
: restrictions(weightedNormalStress.size()) {
auto trivialNonlinearity =
std::make_shared<Friction>(std::make_shared<TrivialFunction>());
: restrictions(), localToGlobal(), zeroFriction() {
auto const gridView = frictionalBoundary.gridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
for (auto it = gridView.template begin<block_size>();
it != gridView.template end<block_size>(); ++it) {
auto const i = vertexMapper.index(*it);
auto const coordinate = it->geometry().corner(0);
if (not frictionalBoundary.containsVertex(i)) {
restrictions[i] = trivialNonlinearity;
if (not frictionalBoundary.containsVertex(i))
continue;
}
auto const fp = std::make_shared<ScalarFriction>(
weights[i], weightedNormalStress[i], frictionInfo(coordinate));
restrictions[i] = std::make_shared<Friction>(fp);
localToGlobal.emplace_back(i);
restrictions.emplace_back(weights[i], weightedNormalStress[i],
frictionInfo(geoToPoint(it->geometry())));
}
assert(restrictions.size() == frictionalBoundary.numVertices());
assert(localToGlobal.size() == frictionalBoundary.numVertices());
}
void updateAlpha(ScalarVector const &alpha) override {
for (size_t i = 0; i < restrictions.size(); ++i)
restrictions[i]->updateAlpha(alpha[i]);
for (size_t j = 0; j < restrictions.size(); ++j)
restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
}
/*
Return a restriction of the outer function to the i'th node.
*/
std::shared_ptr<Friction> restriction(size_t i) const override {
return restrictions[i];
LocalNonlinearity const &restriction(size_t i) const override {
auto const index = indexInSortedRange(localToGlobal, i);
if (index == localToGlobal.size())
return zeroFriction;
return restrictions[index];
}
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
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 @@
template <size_t dimension> class LocalFriction {
public:
virtual ~LocalFriction() {}
using VectorType = Dune::FieldVector<double, dimension>;
using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
explicit LocalFriction(std::shared_ptr<FrictionPotential> func)
: func(func) {}
void virtual updateAlpha(double alpha) = 0;
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 {
return func->evaluate(x.two_norm());
}
void virtual addHessian(VectorType const &x, MatrixType &A) const = 0;
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();
if (xnorm <= 0.0)
return std::numeric_limits<double>::infinity();
return func->regularity(xnorm);
return func_.regularity(xnorm);
}
double coefficientOfFriction(VectorType const &x) const {
return func->coefficientOfFriction(x.two_norm());
double coefficientOfFriction(VectorType const &x) const override {
return func_.coefficientOfFriction(x.two_norm());
}
// directional subdifferential: at u on the line u + t*v
// u and v are assumed to be non-zero
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();
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
D[0] = D[1] = func->differential(xnorm) * (x * v) / xnorm;
D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
}
/** Formula for the derivative:
......@@ -65,14 +85,14 @@ template <size_t dimension> class LocalFriction {
+ \frac {H'(|z|)}{|z|} \operatorname{id}
\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 xnorm = std::sqrt(xnorm2);
if (xnorm2 <= 0.0)
return;
double const H1 = func->differential(xnorm);
double const H2 = func->second_deriv(xnorm);
double const H1 = func_.differential(xnorm);
double const H2 = func_.second_deriv(xnorm);
double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
double const idweight = H1 / xnorm;
......@@ -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();
if (std::isinf(func->regularity(xnorm)))
if (std::isinf(func_.regularity(xnorm)))
return;
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 &,
Dune::Solvers::Interval<double> &dom) const {
Dune::Solvers::Interval<double> &dom) const override {
dom[0] = -std::numeric_limits<double>::max();
dom[1] = std::numeric_limits<double>::max();
}
private:
std::shared_ptr<FrictionPotential> const func;
ScalarFriction func_;
};
#endif
dune/tectonic/minimisation.hh
View file @ 8fe950bc
......@@ -18,7 +18,7 @@ double lineSearch(Functional const &J,
typename Functional::LocalVector const &v,
Bisection const &bisection) {
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;
return bisection.minimize(JRest, 0.0, 0.0, count);
}
......
dune/tectonic/myblockproblem.hh
View file @ 8fe950bc
......@@ -4,7 +4,6 @@
// Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
#include <dune/common/bitsetvector.hh>
#include <dune/common/nullptr.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/fmatrixev.hh>
......@@ -22,16 +21,16 @@
/** \brief Base class for problems where each block can be solved with a
* modified gradient method */
template <class ConvexProblem>
class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
private:
typedef BlockNonlinearGSProblem<ConvexProblem> BNGSP;
typedef BlockNonlinearGSProblem<ConvexProblem> Base;
public:
using typename BNGSP::ConvexProblemType;
using typename BNGSP::LocalMatrixType;
using typename BNGSP::LocalVectorType;
using typename BNGSP::MatrixType;
using typename BNGSP::VectorType;
using typename Base::ConvexProblemType;
using typename Base::LocalMatrixType;
using typename Base::LocalVectorType;
using typename Base::MatrixType;
using typename Base::VectorType;
size_t static const block_size = ConvexProblem::block_size;
size_t static const coarse_block_size = block_size;
......@@ -56,7 +55,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
};
MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
: BNGSP(parset, problem),
: Base(parset, problem),
maxEigenvalues_(problem.f.size()),
localBisection(0.0, 1.0, 0.0, true, 0.0) {
for (size_t i = 0; i < problem.f.size(); ++i) {
......@@ -104,15 +103,10 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
v /= vnorm; // Rescale for numerical stability
MyDirectionalConvexFunction<GlobalFriction<MatrixType, VectorType>> const
psi(computeDirectionalA(problem_.A, v),
computeDirectionalb(problem_.A, problem_.f, u, v), problem_.phi, u,
v);
auto const psi = restrict(problem_.A, problem_.f, u, v, problem_.phi);
Dune::Solvers::Interval<double> D;
psi.subDiff(0, D);
// NOTE: Numerical instability can actually get us here
if (D[1] > 0)
if (D[1] > 0) // NOTE: Numerical instability can actually get us here
return 0;
int bisectionsteps = 0;
......@@ -179,17 +173,15 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
auto const blockEnd = std::end((*col_it)[i]);
for (auto blockIt = std::begin((*col_it)[i]); blockIt != blockEnd;
++blockIt)
if (linearization.truncation[row][i] ||
if (linearization.truncation[row][i] or
linearization.truncation[col][blockIt.index()])
*blockIt = 0.0;
}
}
for (size_t j = 0; j < block_size; ++j)
if (linearization.truncation[row][j])
linearization.b[row][j] = 0.0;
}
for (size_t j = 0; j < block_size; ++j)
outStream << std::setw(9) << linearization.truncation.countmasked(j);
}
......@@ -203,7 +195,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
std::vector<double> maxEigenvalues_;
// problem data
using BNGSP::problem_;
using Base::problem_;
Bisection const localBisection;
......@@ -263,20 +255,17 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
else
ui[j] = 0;
QuadraticEnergy<typename ConvexProblem::NonlinearityType::Friction>
QuadraticEnergy<
typename ConvexProblem::NonlinearityType::LocalNonlinearity>
localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
minimise(localJ, ui, bisection_);
}
}
private:
// problem data
ConvexProblem const &problem;
std::vector<double> maxEigenvalues_;
Bisection const bisection_;
// state data for smoothing procedure used by:
// setIterate, updateIterate, solveLocalProblem
VectorType u;
......
dune/tectonic/mydirectionalconvexfunction.hh
View file @ 8fe950bc
......@@ -7,25 +7,6 @@
#include <dune/fufem/arithmetic.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 {
public:
using Vector = typename Nonlinearity::VectorType;
......@@ -65,4 +46,23 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
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
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 {
using Nonlinearity = NonlinearityTEMPLATE;
using LocalVector = typename Nonlinearity::VectorType;
QuadraticEnergy(double alpha, LocalVector const &b,
std::shared_ptr<Nonlinearity const> phi)
QuadraticEnergy(double alpha, LocalVector const &b, Nonlinearity const &phi)
: alpha(alpha), b(b), phi(phi) {}
double const alpha;
LocalVector const &b;
std::shared_ptr<Nonlinearity const> const phi;
Nonlinearity const &phi;
};
#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
adaptivetimestepper.cc
set(SW_SOURCE_FILES
assemblers.cc
boundary_writer.cc
coupledtimestepper.cc
enumparser.cc
fixedpointiterator.cc
friction_writer.cc
rate.cc
rate/rateupdater.cc
sand-wedge.cc
sand-wedge-data/mygeometry.cc
sand-wedge-data/mygrid.cc
solverfactory.cc
state.cc
hdf5/frictionalboundary-writer.cc
hdf5/iteration-writer.cc
hdf5/patchinfo-writer.cc
hdf5/restart-io.cc
hdf5/surface-writer.cc
hdf5/time-writer.cc
one-body-problem-data/mygeometry.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
)
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)
set(_target sand-wedge-${_dim}D)
add_executable(${_target} ${SOURCE_FILES})
add_dune_pythonlibs_flags(${_target})
add_dune_ug_flags(${_target})
set(_sw_target one-body-problem-${_dim}D)
set(_ugw_target uniform-grid-writer-${_dim}D)
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})
target_link_libraries(${_target} ${Boost_SERIALIZATION_LIBRARY})
target_link_libraries(${_target} ${Boost_SYSTEM_LIBRARY})
add_dune_hdf5_flags(${_sw_target})
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()
src/assemblers.cc
View file @ 8fe950bc
......@@ -7,6 +7,7 @@
#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.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/variablecoefficientviscosityassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
......@@ -14,11 +15,11 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functiontools/p0p1interpolation.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/tectonic/frictionpotential.hh>
#include <dune/tectonic/globalratestatefriction.hh>
#include <dune/tectonic/pointtractionboundaryassembler.hh>
#include "assemblers.hh"
......@@ -111,24 +112,28 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
displacement);
Vector traction;
// Note: We assume v = 0 here so that there is no viscous
// contribution to the stress.
PointTractionBoundaryAssembler<Grid> weightedNormalStressBoundaryAssembler(
Vector traction(cellBasis.size());
NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
youngModulus, poissonRatio, &displacementFunction, 1);
vertexAssembler.assembleBoundaryFunctional(
weightedNormalStressBoundaryAssembler, traction, frictionalBoundary);
std::vector<typename Vector::block_type> normals;
frictionalBoundary.getNormals(normals);
for (size_t i = 0; i < traction.size(); ++i) {
weightedNormalStress[i] = normals[i] * traction[i];
if (weightedNormalStress[i] > 0.0) {
weightedNormalStress[i] = 0.0;
std::cout << "Warning: Manually reducing positive normal stress to zero."
<< std::endl;
}
cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
frictionalBoundary);
auto const nodalTractionAverage =
interpolateP0ToP1(frictionalBoundary, traction);
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>
......@@ -152,13 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
case Config::Truncated:
return std::make_shared<GlobalRateStateFriction<
Matrix, Vector, TruncatedRateState, GridView>>(
frictionalBoundary, gridView, frictionInfo, weights,
weightedNormalStress);
frictionalBoundary, frictionInfo, weights, weightedNormalStress);
case Config::Regularised:
return std::make_shared<GlobalRateStateFriction<
Matrix, Vector, RegularisedRateState, GridView>>(
frictionalBoundary, gridView, frictionInfo, weights,
weightedNormalStress);
frictionalBoundary, frictionInfo, weights, weightedNormalStress);
default:
assert(false);
}
......