diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
index a308c80b63b9cb2b48f0cc096ade7d9cdf464e57..676d1a408726cc0c1a71d3c87010a35182156147 100644
--- a/src/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 0.01 # 2e-3 [m]
+smallestDiameter = 0.5 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
index d58f80cfbe067c712e6d8396d2a37faf5eba7317..77de9404eafc006ed026f47ebfd43b58a4db1fe1 100644
--- a/src/multi-body-problem.cfg
+++ b/src/multi-body-problem.cfg
@@ -11,7 +11,7 @@ vtk.write = true
[problem]
finalTime = 100 # [s] #1000
-bodyCount = 4
+bodyCount = 2
[body]
bulkModulus = 0.5e5 # [Pa]
diff --git a/src/spatial-solving/tnnmg/functional.hh b/src/spatial-solving/tnnmg/functional.hh
index 5a1d3b3736dcb65621342dfdeefeec21f0dfe12c..97dab68876812aefab43d4fe30549bdecf20429e 100644
--- a/src/spatial-solving/tnnmg/functional.hh
+++ b/src/spatial-solving/tnnmg/functional.hh
@@ -16,10 +16,12 @@
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+#include "../../utils/debugutils.hh"
+
template<class M, class V, class N, class L, class U, class R>
class DirectionalRestriction;
-template<class M, class V, class N, class R>
+template<class M, class E, class V, class N, class R>
class ShiftedFunctional;
template <class M, class V, class N, class L, class U, class R>
@@ -33,21 +35,23 @@ class Functional;
* \tparam V nonlinearity type
* \tparam R Range type
*/
-template<class M, class V, class N, class R>
+template<class M, class E, class V, class N, class R>
class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R> {
using Base = Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunctional<M,V,R>;
public:
using Matrix = M;
+ using Eigenvalues = E;
using Vector = V;
using Nonlinearity = N;
using LowerObstacle = Vector;
using UpperObstacle = Vector;
using Range = R;
- ShiftedFunctional(const Matrix& quadraticPart, const Vector& linearPart, const Nonlinearity& phi, const LowerObstacle& lowerObstacle, const UpperObstacle& upperObstacle, const Vector& origin) :
+ ShiftedFunctional(const Matrix& quadraticPart, const Eigenvalues& maxEigenvalues, const Vector& linearPart, const Nonlinearity& phi, const LowerObstacle& lowerObstacle, const UpperObstacle& upperObstacle, const Vector& origin) :
Base(quadraticPart, linearPart, lowerObstacle, upperObstacle, origin),
- phi_(phi)
+ phi_(phi),
+ maxEigenvalues_(maxEigenvalues)
{}
Range operator()(const Vector& v) const
@@ -63,14 +67,13 @@ class ShiftedFunctional : public Dune::TNNMG::ShiftedBoxConstrainedQuadraticFunc
return phi_;
}
- friend auto directionalRestriction(const ShiftedFunctional& f, const Vector& origin, const Vector& direction)
- -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
- {
- return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.originalLinearPart(), f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, direction);
+ const auto& maxEigenvalues() const {
+ return maxEigenvalues_;
}
protected:
const Nonlinearity& phi_;
+ const Eigenvalues& maxEigenvalues_;
};
@@ -167,80 +170,56 @@ class DirectionalRestriction :
* \tparam U Upper obstacle type
* \tparam R Range type
*/
-/*template<class V, class N, class L, class U, class O, class D, class R>
-class FirstOrderModelFunctional {
+template<class N, class V, class R>
+class FirstOrderFunctional {
using Interval = typename Dune::Solvers::Interval<R>;
public:
- using Vector = std::decay_t<V>;
using Nonlinearity = std::decay_t<N>;
- using LowerObstacle = std::decay_t<L>;
- using UpperObstacle = std::decay_t<U>;
+ using Vector = V;
+ using LowerObstacle = V;
+ using UpperObstacle = V;
using Range = R;
- using field_type = typename Vector::field_type;
-
+ using field_type = typename V::field_type;
public:
- template <class VV1, class NN, class LL, class UU, class VV2>
- FirstOrderModelFunctional(
+ template <class NN, class LL, class UU, class OO, class DD>
+ FirstOrderFunctional(
const Range& maxEigenvalue,
- VV1&& linearTerm,
+ const Range& linearPart,
NN&& phi,
LL&& lower,
UU&& upper,
- VV2&& origin) :
- maxEigenvalue_(maxEigenvalue),
- linearTerm_(std::forward<VV1>(linearTerm)),
+ OO&& origin,
+ DD&& direction) :
+ quadraticPart_(maxEigenvalue),
+ linearPart_(linearPart),
lower_(std::forward<LL>(lower)),
upper_(std::forward<UU>(upper)),
- origin_(std::forward<VV2>(origin)),
- //direction_(linearTerm_),
+ origin_(std::forward<OO>(origin)),
+ direction_(std::forward<DD>(direction)),
phi_(std::forward<NN>(phi)) {
- auto& origin = origin_.get();
- auto& direction = direction_.get();
-
- auto v = ri;
- double const vnorm = v.two_norm();
- if (vnorm > 1.0)
- v /= vnorm;
-
// set defect obstacles
defectLower_ = -std::numeric_limits<field_type>::max();
defectUpper_ = std::numeric_limits<field_type>::max();
- Dune::TNNMG::directionalObstacles(origin, direction, lower_.get(), upper_.get(), defectLower_, defectUpper_);
+ Dune::TNNMG::directionalObstacles(origin_.get(), direction_.get(), lower_.get(), upper_.get(), defectLower_, defectUpper_);
// set domain
- phi_.get().directionalDomain(origin, direction, domain_);
+ phi_.get().directionalDomain(origin_.get(), direction_.get(), domain_);
- if (domain_[0] < defectLower_) {
+ /*if (domain_[0] < defectLower_) {
domain_[0] = defectLower_;
}
if (domain_[1] > defectUpper_) {
domain_[1] = defectUpper_;
- }
+ }*/
}
- template <class BitVector>
- auto setIgnore(const BitVector& ignore) const {
- Vector direction = direction_.get();
- Vector origin = origin_.get();
- // Dune::Solvers::resizeInitializeZero(direction, linearPart());
- for (size_t j = 0; j < ignore.size(); ++j)
- if (ignore[j])
- direction[j] = 0; // makes sure result remains in subspace after correction
- else
- origin[j] = 0; // shift affine offset
-
- // set quadratic and linear parts of functional
- quadraticPart_ = maxEigenvalue_ * direction.two_norm2();
- linearPart_ = linearTerm_.get() * direction;
- }
-
Range operator()(const Vector& v) const
{
- DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderModelFunctional not implemented");
+ DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderFunctional not implemented");
}
auto subDifferential(double x) const {
@@ -261,49 +240,53 @@ class FirstOrderModelFunctional {
return domain_;
}
- const Vector& origin() const {
+ const auto& origin() const {
return origin_.get();
}
- const Vector& direction() const {
+ const auto& direction() const {
return direction_.get();
}
- const field_type& defectLower() const {
+ const auto& lowerObstacle() const {
return defectLower_;
}
- const field_type& defectUpper() const {
+ const auto& upperObstacle() const {
return defectUpper_;
}
+ const auto& quadraticPart() const {
+ return quadraticPart_;
+ }
+
+ const auto& linearPart() const {
+ return linearPart_;
+ }
private:
- Dune::Solvers::ConstCopyOrReference<V> linearTerm_;
- Dune::Solvers::ConstCopyOrReference<L> lower_;
- Dune::Solvers::ConstCopyOrReference<U> upper_;
+ const Range quadraticPart_;
+ const Range linearPart_;
+
+ Dune::Solvers::ConstCopyOrReference<LowerObstacle> lower_;
+ Dune::Solvers::ConstCopyOrReference<UpperObstacle> upper_;
- Dune::Solvers::ConstCopyOrReference<O> origin_;
- Dune::Solvers::ConstCopyOrReference<D> direction_;
+ Dune::Solvers::ConstCopyOrReference<Vector> origin_;
+ Dune::Solvers::ConstCopyOrReference<Vector> direction_;
Dune::Solvers::ConstCopyOrReference<N> phi_;
Interval domain_;
- Range maxEigenvalue_;
-
- field_type quadraticPart_;
- field_type linearPart_;
-
field_type defectLower_;
field_type defectUpper_;
-};*/
+};
// \ToDo This should be an inline friend of ShiftedBoxConstrainedQuadraticFunctional
// but gcc-4.9.2 shipped with debian jessie does not accept
// inline friends with auto return type due to bug-59766.
// Notice, that this is fixed in gcc-4.9.3.
-template<class Index, class M, class V, class Nonlinearity, class R>
-auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, const Index& i)
+template<class Index, class M, class E, class V, class Nonlinearity, class R>
+auto coordinateRestriction(const ShiftedFunctional<M, E, V, Nonlinearity, R>& f, const Index& i)
{
using Range = R;
using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
@@ -316,6 +299,8 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
const LocalMatrix* Aii_p = nullptr;
+ //print(f.origin(), "f.origin: ");
+
LocalVector ri = f.originalLinearPart()[i];
const auto& Ai = f.quadraticPart()[i];
sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
@@ -326,35 +311,9 @@ auto coordinateRestriction(const ShiftedFunctional<M, V, Nonlinearity, R>& f, co
Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
});
- // set maxEigenvalue from matrix
- LocalVector eigenvalues;
- Dune::FMatrixHelp::eigenValues(*Aii_p, eigenvalues);
- auto maxEigenvalue = *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
-
- Dune::MatrixVector::addProduct(ri, maxEigenvalue, f.origin()[i]);
-
- LocalMatrix Aii = 0;
- for (size_t d=0; d<Aii.N(); d++) {
- Aii[d][d] = maxEigenvalue;
- }
-
- LocalLowerObstacle dli = f.originalLowerObstacle()[i];
- LocalUpperObstacle dui = f.originalUpperObstacle()[i];
- dli -= f.origin()[i];
- dui -= f.origin()[i];
-
auto&& phii = f.phi().restriction(i);
- /*std::cout << "maxEigenvalue matrix: " << std::endl;
- for (size_t i=0; i<Aii.N(); i++) {
- for (size_t j=0; j<Aii.M(); j++) {
- std::cout << Aii[i][j] << " ";
- }
- std::cout << std::endl;
- }*/
-
- //return FirstOrderModelFunctional<LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, LocalVector, Range>(maxEigenvalue, std::move(ri), std::move(phii), std::move(dli), std::move(dui), std::move(f.origin()[i]));
- return ShiftedFunctional<LocalMatrix, LocalVector, decltype(phii), Range>(std::move(Aii), std::move(ri), std::move(phii), std::move(dli), std::move(dui), f.origin()[i]);
+ return ShiftedFunctional<LocalMatrix, Range, LocalVector, decltype(phii), Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), std::move(phii), f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
}
@@ -382,8 +341,11 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
using LowerObstacle = typename Base::LowerObstacle;
using UpperObstacle = typename Base::UpperObstacle;
+ using Eigenvalues = std::vector<Range>;
+
private:
Dune::Solvers::ConstCopyOrReference<N> phi_;
+ Eigenvalues maxEigenvalues_;
public:
@@ -395,13 +357,25 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
LL&& lower,
UU&& upper) :
Base(std::forward<MM>(matrix), std::forward<VV>(linearPart), std::forward<LL>(lower), std::forward<UU>(upper)),
- phi_(std::forward<NN>(phi))
- {}
+ phi_(std::forward<NN>(phi)),
+ maxEigenvalues_(this->linearPart().size())
+ {
+ for (size_t i=0; i<this->quadraticPart().N(); ++i) {
+ typename Vector::block_type eigenvalues;
+ Dune::FMatrixHelp::eigenValues(this->quadraticPart()[i][i], eigenvalues);
+ maxEigenvalues_[i] =
+ *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+ }
+ }
const Nonlinearity& phi() const {
return phi_.get();
}
+ const auto& maxEigenvalues() const {
+ return maxEigenvalues_;
+ }
+
Range operator()(const Vector& v) const
{
if (Dune::TNNMG::checkInsideBox(v, this->lower_.get(), this->upper_.get()))
@@ -416,9 +390,9 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
}
friend auto shift(const Functional& f, const Vector& origin)
- -> ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>
+ -> ShiftedFunctional<Matrix, Eigenvalues, Vector, Nonlinearity, Range>
{
- return ShiftedFunctional<Matrix, Vector, Nonlinearity, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+ return ShiftedFunctional<Matrix, Eigenvalues, Vector, Nonlinearity, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
}
};
diff --git a/src/spatial-solving/tnnmg/localbisectionsolver.hh b/src/spatial-solving/tnnmg/localbisectionsolver.hh
index 33f69e956c5c0e6b9102164932d4a440b08c611b..b24eb8ee6f84da9652d11924ea1a44056a1836c8 100644
--- a/src/spatial-solving/tnnmg/localbisectionsolver.hh
+++ b/src/spatial-solving/tnnmg/localbisectionsolver.hh
@@ -3,8 +3,13 @@
#ifndef DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
#define DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "functional.hh"
+
#include "../../utils/debugutils.hh"
/**
@@ -21,26 +26,33 @@ class LocalBisectionSolver
return;
}
- Vector origin = f.origin();
- Vector direction = f.originalLinearPart();
+ auto maxEigenvalue = f.maxEigenvalues();
+
+ auto origin = f.origin();
+ auto linearPart = f.originalLinearPart();
+ Dune::MatrixVector::addProduct(linearPart, maxEigenvalue, origin);
for (size_t j = 0; j < ignore.size(); ++j)
if (ignore[j])
- direction[j] = 0; // makes sure result remains in subspace after correction
+ linearPart[j] = 0; // makes sure result remains in subspace after correction
else
origin[j] = 0; // shift affine offset
- double const directionNorm = direction.two_norm();
- if (directionNorm <= 0.0)
+ double const linearPartNorm = linearPart.two_norm();
+ if (linearPartNorm <= 0.0)
return;
- direction /= directionNorm;
+ linearPart /= linearPartNorm;
- auto&& directionalF = directionalRestriction(f, origin, direction);
+ using Nonlinearity = std::decay_t<decltype(f.phi())>;
+ using Range = std::decay_t<decltype(f.maxEigenvalues())>;
+ using FirstOrderFunctional = FirstOrderFunctional<Nonlinearity, Vector, Range>;
+
+ FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, linearPart);
// scalar obstacle solver without nonlinearity
typename Functional::Range alpha;
- Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
- obstacleSolver(alpha, directionalF, false);
+ /*Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+ obstacleSolver(alpha, firstOrderFunctional, false);*/
//direction *= alpha;
@@ -59,23 +71,31 @@ class LocalBisectionSolver
auto D = directionalF.subDifferential(0);
std::cout << "subDiff at x: " << D[0] << " " << D[1] << std::endl;*/
- const auto& domain = directionalF.domain();
+ const auto& domain = firstOrderFunctional.domain();
if (std::abs(domain[0]-domain[1])>safety) {
int bisectionsteps = 0;
- const Bisection globalBisection;
+ const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);
+
+ alpha = globalBisection.minimize(firstOrderFunctional, 0.0, 0.0, bisectionsteps);
+
+ if (alpha < firstOrderFunctional.lowerObstacle() ) {
+ alpha = firstOrderFunctional.lowerObstacle();
+ }
- alpha = globalBisection.minimize(directionalF, alpha, 0.0, bisectionsteps);
+ if (alpha > firstOrderFunctional.upperObstacle() ) {
+ alpha = firstOrderFunctional.upperObstacle();
+ }
}
- direction *= alpha;
+ linearPart *= alpha;
#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
x = origin;
- x += direction;
+ x += linearPart;
#else
//x = origin;
- x -= f.origin();
- x += direction;
+ //x -= f.origin();
+ x = linearPart;
#endif
}
};