Skip to content
Snippets Groups Projects

Compare revisions

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

Source

Select target project
No results found

Target

Select target project
  • podlesny/dune-tectonic
  • agnumpde/dune-tectonic
2 results
Show changes
Commits on Source (89)
Hide whitespace changes
Inline Side-by-side
Showing
with 260 additions and 717 deletions
CMakeLists.txt
View file @ 8fe950bc
......@@ -18,6 +18,8 @@ include(DuneMacros)
# 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
#Libs: -L${libdir} -ldunetectonic
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
......@@ -2,9 +2,7 @@
# Dune module information file #
################################
#Name of the module
Module: dune-tectonic
Version: 2.1-1
Maintainer: pipping@mi.fu-berlin.de
#depending on
Depends: dune-common dune-fufem dune-tnnmg
Version: 2.5-1
Maintainer: elias.pipping@fu-berlin.de
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
......@@ -29,35 +29,36 @@ class FrictionPotential {
class TruncatedRateState : public FrictionPotential {
public:
TruncatedRateState(double coefficient, double _normalStress, FrictionData _fd)
: fd(_fd), weight(coefficient), normalStress(_normalStress) {}
TruncatedRateState(double _weight, double _weightedNormalStress,
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 {
return weight * (fd.C - normalStress * coefficientOfFriction(V));
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 weight * (-normalStress) * (fd.a / V);
return -weightedNormalStress * (fd.a / V);
}
double regularity(double V) const {
double regularity(double V) const override {
if (std::abs(V - Vmin) < 1e-14) // TODO
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);
}
......@@ -65,31 +66,33 @@ class TruncatedRateState : public FrictionPotential {
private:
FrictionData const fd;
double const weight;
double const normalStress;
double const weightedNormalStress;
double Vmin;
};
class RegularisedRateState : public FrictionPotential {
public:
RegularisedRateState(double coefficient, double _normalStress,
RegularisedRateState(double _weight, double _weightedNormalStress,
FrictionData _fd)
: fd(_fd), weight(coefficient), normalStress(_normalStress) {}
: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
double coefficientOfFriction(double V) const {
double coefficientOfFriction(double V) const override {
return fd.a * std::asinh(0.5 * V / Vmin);
}
double differential(double V) const {
return weight * (fd.C - normalStress * coefficientOfFriction(V));
double differential(double V) const override {
return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
}
double second_deriv(double V) const {
return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V);
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);
}
......@@ -97,22 +100,22 @@ class RegularisedRateState : public FrictionPotential {
private:
FrictionData const fd;
double const weight;
double const normalStress;
double const weightedNormalStress;
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) {
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,58 +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,
// Note: passing the following two makes no sense
ScalarVector const &weights,
ScalarVector const &normalStress)
: restrictions(normalStress.size()) {
auto trivialNonlinearity =
std::make_shared<Friction>(std::make_shared<TrivialFunction>());
ScalarVector const &weightedNormalStress)
: 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], normalStress[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/polyhedrondistance.hh deleted 100644 → 0
View file @ 4bea0810
#ifndef DUNE_TECTONIC_POLYHEDRONDISTANCE_HH
#define DUNE_TECTONIC_POLYHEDRONDISTANCE_HH
// Based on the closest point projection from dune-contact
#include <dune/common/dynmatrix.hh>
#include <dune/solvers/common/arithmetic.hh>
template <class Coordinate> struct ConvexPolyhedron {
std::vector<Coordinate> vertices;
template <class DynamicVector>
Coordinate barycentricToGlobal(DynamicVector const &b) const {
assert(b.size() == vertices.size());
Coordinate r(0);
for (size_t i = 0; i < b.size(); i++)
Arithmetic::addProduct(r, b[i], vertices[i]);
return r;
}
template <class DynamicVector>
void sanityCheck(DynamicVector const &b, double tol) const {
double sum = 0;
for (size_t i = 0; i < b.size(); i++) {
if (b[i] < -tol)
DUNE_THROW(Dune::Exception, "No barycentric coords " << b[1] << " nr. "
<< i);
sum += b[i];
}
if (sum > 1 + tol)
DUNE_THROW(Dune::Exception, "No barycentric coords, sum: " << sum);
}
};
template <class Coordinate, class Matrix>
void populateMatrix(ConvexPolyhedron<Coordinate> const &segment,
Matrix &matrix) {
size_t const nCorners = segment.vertices.size();
for (size_t i = 0; i < nCorners; i++)
for (size_t j = 0; j <= i; j++)
matrix[i][j] = segment.vertices[i] * segment.vertices[j];
// make symmetric
for (size_t i = 0; i < nCorners; i++)
for (size_t j = i + 1; j < nCorners; j++)
matrix[i][j] = matrix[j][i];
}
/**
* The functional to be minimized is
*
* 0.5*norm(\sum_i lambda_i corner_i - target)^2
*
* where lambda_i are barycentric coordinates, i.e.
*
* 0 \le lambda_i \le 1 and \sum_i lambda_i=1
*
* The resulting quadratic minimization problem is given by
*
* 0.5 lambda^T*A*lambda - l^T*lambda
*
* with A_ij = (corner_i,corner_j) and l_i = (corner_i, target)
*/
template <class Coordinate>
void approximate(Coordinate const &target, Dune::DynamicMatrix<double> const &A,
Dune::DynamicVector<double> const &l,
ConvexPolyhedron<Coordinate> const &segment,
Dune::DynamicVector<double> &sol) {
size_t nCorners = segment.vertices.size();
// compute the residual
Dune::DynamicVector<double> rhs = l;
// compute rhs -= A*sol_
for (size_t k = 0; k < nCorners; k++)
for (size_t j = 0; j < nCorners; j++)
rhs[k] -= A[k][j] * sol[j];
// use the edge vectors as search directions
double alpha = 0.0;
for (size_t j1 = 0; j1 < nCorners - 1; ++j1) {
for (size_t j2 = j1 + 1; j2 < nCorners; ++j2) {
// compute matrix entry and rhs for edge direction
double a = A[j1][j1] - A[j1][j2] - A[j2][j1] + A[j2][j2];
double b = rhs[j1] - rhs[j2];
// compute minimizer for correction
if (a > 0) {
alpha = b / a;
double const largestAlpha = sol[j2];
double const smallestAlpha = -sol[j1];
// project alpha such that we stay positive
if (alpha > largestAlpha)
alpha = largestAlpha;
else if (alpha < smallestAlpha)
alpha = smallestAlpha;
} else {
// if the quadratic term vanishes, we have a linear function, which
// attains its extrema on the boundary
double sum = sol[j1] + sol[j2];
double lValue = 0.5 * A[j1][j1] * sum * sum - rhs[j1] * sum;
double uValue = 0.5 * A[j2][j2] * sum * sum - rhs[j2] * sum;
alpha = lValue < uValue ? sol[j2] : -sol[j1];
}
// apply correction
sol[j1] += alpha;
sol[j2] -= alpha;
// update the local residual for corrections
for (size_t p = 0; p < nCorners; p++)
rhs[p] -= (A[p][j1] - A[p][j2]) * alpha;
}
}
}
// Point-to-Polyhedron
template <class Coordinate>
double distance(const Coordinate &target,
const ConvexPolyhedron<Coordinate> &segment,
double correctionTolerance) {
size_t nCorners = segment.vertices.size();
double const barycentricTolerance = 1e-14;
Dune::DynamicMatrix<double> A(nCorners, nCorners);
populateMatrix(segment, A);
Dune::DynamicVector<double> l(nCorners);
for (size_t i = 0; i < nCorners; i++)
l[i] = target * segment.vertices[i];
// choose a feasible initial solution
Dune::DynamicVector<double> sol(nCorners);
for (size_t i = 0; i < nCorners; i++)
sol[i] = 1.0 / nCorners;
Coordinate oldCoordinates = segment.barycentricToGlobal(sol);
Coordinate newCoordinates;
size_t maxIterations = 1000;
size_t iterations;
for (iterations = 0; iterations < maxIterations; ++iterations) {
approximate(target, A, l, segment, sol);
newCoordinates = segment.barycentricToGlobal(sol);
if ((oldCoordinates - newCoordinates).two_norm() < correctionTolerance)
break;
oldCoordinates = newCoordinates;
}
assert(iterations < maxIterations);
segment.sanityCheck(sol, barycentricTolerance);
return (target - newCoordinates).two_norm();
}
// Polyhedron-to-Polyhedron
template <class Coordinate>
double distance(const ConvexPolyhedron<Coordinate> &s1,
const ConvexPolyhedron<Coordinate> &s2,
double valueCorrectionTolerance) {
size_t nCorners1 = s1.vertices.size();
size_t nCorners2 = s2.vertices.size();
double const barycentricTolerance = 1e-14;
// choose feasible initial solutions
Dune::DynamicVector<double> sol1(nCorners1);
for (size_t i = 0; i < s1.vertices.size(); i++)
sol1[i] = 1.0 / s1.vertices.size();
auto c1 = s1.barycentricToGlobal(sol1);
Dune::DynamicVector<double> sol2(nCorners2);
for (size_t i = 0; i < nCorners2; i++)
sol2[i] = 1.0 / nCorners2;
auto c2 = s2.barycentricToGlobal(sol2);
double oldDistance = (c2 - c1).two_norm();
Dune::DynamicMatrix<double> A1(nCorners1, nCorners1);
populateMatrix(s1, A1);
Dune::DynamicMatrix<double> A2(nCorners2, nCorners2);
populateMatrix(s2, A2);
size_t maxIterations = 1000;
size_t iterations;
for (iterations = 0; iterations < maxIterations; ++iterations) {
Dune::DynamicVector<double> l2(nCorners2);
for (size_t i = 0; i < nCorners2; i++)
l2[i] = c1 * s2.vertices[i];
approximate(c1, A2, l2, s2, sol2);
c2 = s2.barycentricToGlobal(sol2);
Dune::DynamicVector<double> l1(nCorners1);
for (size_t i = 0; i < nCorners1; i++)
l1[i] = c2 * s1.vertices[i];
approximate(c2, A1, l1, s1, sol1);
c1 = s1.barycentricToGlobal(sol1);
double const newDistance = (c2 - c1).two_norm();
assert(newDistance - oldDistance <= 1e-12); // debugging
if (oldDistance - newDistance <= valueCorrectionTolerance) {
s1.sanityCheck(sol1, barycentricTolerance);
s2.sanityCheck(sol2, barycentricTolerance);
return newDistance;
}
oldDistance = newDistance;
}
assert(false);
}
#endif
dune/tectonic/quadraticenergy.hh
View file @ 8fe950bc
......@@ -8,12 +8,11 @@ template <class NonlinearityTEMPLATE> class QuadraticEnergy {
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 @ 4bea0810
set(TESTS test-polyhedral-minimisation)
foreach(_test ${TESTS})
add_executable(${_test} EXCLUDE_FROM_ALL ${_test}.cc)
target_link_dune_default_libraries(${_test})
add_test(${_test} ${_test})
endforeach()
add_directory_test_target(_test_target)
add_dependencies(${_test_target} ${TESTS})
dune/tectonic/test/test-polyhedral-minimisation.cc deleted 100644 → 0
View file @ 4bea0810
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/fvector.hh>
#include <dune/tectonic/polyhedrondistance.hh>
int main() {
using LocalVector = Dune::FieldVector<double, 2>;
auto const test =
[&](ConvexPolyhedron<LocalVector> const &p1,
ConvexPolyhedron<LocalVector> const &p2, double analyticalDistance) {
LocalVector target;
{
double const error =
std::abs(analyticalDistance - distance(p1, p2, 1e-12));
std::cout << "error: " << error << std::endl;
assert(error < 1e-12);
}
};
{
/*
* Calculate the distance between two triangles, where it is attained at a
* face-vertex pair
*
* O
* |\
* | \
* O O--O
* |\
* | \
* O--O
*/
double const analyticalDistance = std::sqrt(2.0);
LocalVector tmp;
ConvexPolyhedron<LocalVector> bs1;
bs1.vertices.resize(3);
tmp[0] = 0;
tmp[1] = 0;
bs1.vertices[0] = tmp;
tmp[0] = 0;
tmp[1] = 2;
bs1.vertices[1] = tmp;
tmp[0] = 2;
tmp[1] = 0;
bs1.vertices[2] = tmp;
ConvexPolyhedron<LocalVector> bs2;
bs2.vertices.resize(3);
tmp[0] = 2;
tmp[1] = 2;
bs2.vertices[0] = tmp;
tmp[0] = 2;
tmp[1] = 4;
bs2.vertices[1] = tmp;
tmp[0] = 4;
tmp[1] = 2;
bs2.vertices[2] = tmp;
test(bs1, bs2, analyticalDistance);
}
{
/*
* Calculate the distance between two triangles, where it is
* attained in a face-face pair
*
* O--O
* \ |
* \|
* O O
* |\
* | \
* O--O
*/
double const analyticalDistance = 2.0 * std::sqrt(2.0);
LocalVector tmp;
ConvexPolyhedron<LocalVector> bs1;
bs1.vertices.resize(3);
tmp[0] = 0;
tmp[1] = 0;
bs1.vertices[0] = tmp;
tmp[0] = 0;
tmp[1] = 2;
bs1.vertices[1] = tmp;
tmp[0] = 2;
tmp[1] = 0;
bs1.vertices[2] = tmp;
ConvexPolyhedron<LocalVector> bs2;
bs2.vertices.resize(3);
tmp[0] = 4;
tmp[1] = 4;
bs2.vertices[0] = tmp;
tmp[0] = 2;
tmp[1] = 4;
bs2.vertices[1] = tmp;
tmp[0] = 4;
tmp[1] = 2;
bs2.vertices[2] = tmp;
test(bs1, bs2, analyticalDistance);
}
{
/*
* Calculate the distance between two triangles, where it is
* attained in a vertex-vertex pair
*
* O
* |\
* | \
* O--O O--O
* \ |
* \|
* O
*/
double analyticalDistance = 2.0;
LocalVector tmp;
ConvexPolyhedron<LocalVector> bs1;
bs1.vertices.resize(3);
tmp[0] = 2;
tmp[1] = 2;
bs1.vertices[0] = tmp;
tmp[0] = 0;
tmp[1] = 2;
bs1.vertices[1] = tmp;
tmp[0] = 2;
tmp[1] = 0;
bs1.vertices[2] = tmp;
ConvexPolyhedron<LocalVector> bs2;
bs2.vertices.resize(3);
tmp[0] = 4;
tmp[1] = 2;
bs2.vertices[0] = tmp;
tmp[0] = 4;
tmp[1] = 4;
bs2.vertices[1] = tmp;
tmp[0] = 6;
tmp[1] = 2;
bs2.vertices[2] = tmp;
test(bs1, bs2, analyticalDistance);
}
}
src/CMakeLists.txt
View file @ 8fe950bc
set(SOURCE_FILES
set(SW_SOURCE_FILES
assemblers.cc
boundary_writer.cc
coupledtimestepper.cc
enumparser.cc
fixedpointiterator.cc
friction_writer.cc
sand-wedge.cc
sand-wedge-data/mygeometry.cc
sand-wedge-data/mygrid.cc
solverfactory.cc
state.cc
timestepping.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)
include_directories(${Boost_INCLUDE_DIR})
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_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/adaptivetimestepper.hh deleted 100644 → 0
View file @ 4bea0810
#include "coupledtimestepper.hh"
template <typename T1, typename T2>
std::pair<T1, T2> clonePair(std::pair<T1, T2> in) {
return { in.first->clone(), in.second->clone() };
}
template <class Factory, class UpdaterPair, class ErrorNorm>
class AdaptiveTimeStepper {
using StateUpdater = typename UpdaterPair::first_type::element_type;
using VelocityUpdater = typename UpdaterPair::second_type::element_type;
using Vector = typename Factory::Vector;
using ConvexProblem = typename Factory::ConvexProblem;
using Nonlinearity = typename ConvexProblem::NonlinearityType;
using MyCoupledTimeStepper =
CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>;
public:
AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, UpdaterPair &current,
std::function<void(double, Vector &)> externalForces,
ErrorNorm const &errorNorm,
std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine)
: finalTime_(parset.get<double>("problem.finalTime")),
relativeTime_(0.0),
relativeTau_(1e-6), // FIXME (not really important, though)
factory_(factory),
parset_(parset),
globalFriction_(globalFriction),
current_(current),
R1_(clonePair(current_)),
externalForces_(externalForces),
mustRefine_(mustRefine),
iterationWriter_("iterations", std::fstream::out),
errorNorm_(errorNorm) {
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, R1_.first, R1_.second,
errorNorm, externalForces_);
stepAndReport("R1", coupledTimeStepper, relativeTime_, relativeTau_);
iterationWriter_ << std::endl;
}
bool reachedEnd() { return relativeTime_ >= 1.0; }
bool coarsen() {
bool didCoarsen = false;
while (relativeTime_ + relativeTau_ < 1.0) {
R2_ = clonePair(R1_);
{
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, R2_.first,
R2_.second, errorNorm_, externalForces_);
stepAndReport("R2", coupledTimeStepper, relativeTime_ + relativeTau_,
relativeTau_);
}
UpdaterPair C = clonePair(current_);
{
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, C.first, C.second,
errorNorm_, externalForces_);
stepAndReport("C", coupledTimeStepper, relativeTime_,
2.0 * relativeTau_);
}
if (!mustRefine_(C, R2_)) {
R2_ = { nullptr, nullptr };
R1_ = C;
relativeTau_ *= 2.0;
didCoarsen = true;
} else {
break;
}
}
return didCoarsen;
}
void refine() {
while (true) {
UpdaterPair F2 = clonePair(current_);
UpdaterPair F1;
{
MyCoupledTimeStepper coupledTimeStepper(
finalTime_, factory_, parset_, globalFriction_, F2.first, F2.second,
errorNorm_, externalForces_);
stepAndReport("F1", coupledTimeStepper, relativeTime_,
relativeTau_ / 2.0);
F1 = clonePair(F2);
stepAndReport("F2", coupledTimeStepper,
relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
}
if (!mustRefine_(R1_, F2)) {
break;
} else {
R1_ = F1;
R2_ = F2;
relativeTau_ /= 2.0;
}
}
}
void advance() {
if (!coarsen())
refine();
iterationWriter_ << std::endl;
current_ = R1_;
R1_ = R2_;
relativeTime_ += relativeTau_;
}
double getRelativeTime() { return relativeTime_; }
double getRelativeTau() { return relativeTau_; }
private:
void stepAndReport(std::string type, MyCoupledTimeStepper &stepper,
double rTime, double rTau) {
iterationWriter_ << type << " " << stepper.step(rTime, rTau) << " "
<< std::flush;
}
double finalTime_;
double relativeTime_;
double relativeTau_;
Factory &factory_;
Dune::ParameterTree const &parset_;
std::shared_ptr<Nonlinearity> globalFriction_;
UpdaterPair &current_;
UpdaterPair R1_;
UpdaterPair R2_;
std::function<void(double, Vector &)> externalForces_;
std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine_;
std::fstream iterationWriter_;
ErrorNorm const &errorNorm_;
};
src/assemblers.cc
View file @ 8fe950bc
......@@ -7,14 +7,15 @@
#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>
#include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/computestress.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>
......@@ -33,7 +34,7 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
BoundaryPatch<GridView> const &frictionalBoundary,
ScalarMatrix &frictionalBoundaryMass) {
ScalarMatrix &frictionalBoundaryMass) const {
BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
frictionalBoundaryMassAssembler(frictionalBoundary);
......@@ -45,7 +46,7 @@ template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleMass(
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
densityFunction,
Matrix &M) {
Matrix &M) const {
// NOTE: We treat the weight as a constant function
QuadratureRuleKey quadKey(dimension, 0);
......@@ -58,7 +59,7 @@ void MyAssembler<GridView, dimension>::assembleMass(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
Matrix &A) {
Matrix &A) const {
StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
localStiffness(E, nu);
vertexAssembler.assembleOperator(localStiffness, A);
......@@ -68,7 +69,7 @@ template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleViscosity(
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
Matrix &C) {
Matrix &C) const {
// NOTE: We treat the weights as constant functions
QuadratureRuleKey shearViscosityKey(dimension, 0);
QuadratureRuleKey bulkViscosityKey(dimension, 0);
......@@ -83,7 +84,7 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleBodyForce(
Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
Vector &f) {
Vector &f) const {
L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
gravityFunctionalAssembler(gravityField);
vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
......@@ -92,7 +93,8 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNeumann(
BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
Dune::VirtualFunction<double, double> const &neumann, double relativeTime) {
Dune::VirtualFunction<double, double> const &neumann,
double relativeTime) const {
LocalVector localNeumann(0);
neumann.evaluate(relativeTime, localNeumann[0]);
NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
......@@ -103,25 +105,35 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
}
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNormalStress(
void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
BoundaryPatch<GridView> const &frictionalBoundary,
ScalarVector &normalStress, double youngModulus, double poissonRatio,
Vector const &displacement) {
Vector traction;
Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!)
(frictionalBoundary, displacement, traction, youngModulus, poissonRatio);
std::vector<typename Vector::block_type> normals;
frictionalBoundary.getNormals(normals);
for (size_t i = 0; i < traction.size(); ++i) {
normalStress[i] = normals[i] * traction[i];
if (normalStress[i] > 0.0) {
normalStress[i] = 0.0;
std::cout << "Warning: Manually reducing positive normal stress to zero."
<< std::endl;
}
ScalarVector &weightedNormalStress, double youngModulus,
double poissonRatio, Vector const &displacement) const {
BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
displacement);
Vector traction(cellBasis.size());
NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
youngModulus, poissonRatio, &displacementFunction, 1);
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>
......@@ -129,9 +141,9 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
Config::FrictionModel frictionModel,
BoundaryPatch<GridView> const &frictionalBoundary,
GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &normalStress)
ScalarVector const &weightedNormalStress) const
-> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
// Lump negative normal stress (kludge)
// Lumping of the nonlinearity
ScalarVector weights;
{
NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
......@@ -145,11 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
case Config::Truncated:
return std::make_shared<GlobalRateStateFriction<
Matrix, Vector, TruncatedRateState, GridView>>(
frictionalBoundary, gridView, frictionInfo, weights, normalStress);
frictionalBoundary, frictionInfo, weights, weightedNormalStress);
case Config::Regularised:
return std::make_shared<GlobalRateStateFriction<
Matrix, Vector, RegularisedRateState, GridView>>(
frictionalBoundary, gridView, frictionInfo, weights, normalStress);
frictionalBoundary, frictionInfo, weights, weightedNormalStress);
default:
assert(false);
}
......@@ -158,7 +170,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleVonMisesStress(
double youngModulus, double poissonRatio, Vector const &u,
ScalarVector &stress) {
ScalarVector &stress) const {
auto const gridDisplacement =
std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
vertexBasis, u);
......