Skip to content
Snippets Groups Projects
Commit 569cf70d authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

LocalNonlinearity -> LocalFriction

Also kill the kink logic and thus minimisation2
parent 7d64eb7d
No related branches found
No related tags found
No related merge requests found
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
#include <dune/fufem/interval.hh> #include <dune/fufem/interval.hh>
#include <dune/fufem/arithmetic.hh> #include <dune/fufem/arithmetic.hh>
#include "localnonlinearity.hh" #include "localfriction.hh"
namespace Dune { namespace Dune {
template <int dim> class EllipticEnergy { template <int dim> class EllipticEnergy {
...@@ -16,14 +16,12 @@ template <int dim> class EllipticEnergy { ...@@ -16,14 +16,12 @@ template <int dim> class EllipticEnergy {
using SmallVector = FieldVector<double, dim>; using SmallVector = FieldVector<double, dim>;
using SmallMatrix = FieldMatrix<double, dim, dim>; using SmallMatrix = FieldMatrix<double, dim, dim>;
using NonlinearityType = LocalNonlinearity<dim>; using NonlinearityType = LocalFriction<dim>;
EllipticEnergy(SmallMatrix const &A, SmallVector const &b, EllipticEnergy(SmallMatrix const &A, SmallVector const &b,
shared_ptr<NonlinearityType const> phi, int ignore = dim) shared_ptr<NonlinearityType const> phi, int ignore = dim)
: A(A), b(b), phi(phi), ignore(ignore) {} : A(A), b(b), phi(phi), ignore(ignore) {}
std::vector<double> const &get_kinks() const { return phi->get_kinks(); }
double operator()(SmallVector const &v) const { double operator()(SmallVector const &v) const {
SmallVector y(0); SmallVector y(0);
Arithmetic::addProduct(y, 0.5, A, v); Arithmetic::addProduct(y, 0.5, A, v);
......
...@@ -13,17 +13,7 @@ ...@@ -13,17 +13,7 @@
namespace Dune { namespace Dune {
class FrictionPotentialWrapper { class FrictionPotentialWrapper {
protected:
FrictionPotentialWrapper(std::vector<double> const &_kinks) : kinks(_kinks) {}
private:
std::vector<double> const kinks;
public: public:
std::vector<double> const &get_kinks() const { return kinks; }
FrictionPotentialWrapper() : kinks() {}
virtual ~FrictionPotentialWrapper() {} virtual ~FrictionPotentialWrapper() {}
double virtual leftDifferential(double s) const = 0; double virtual leftDifferential(double s) const = 0;
...@@ -112,8 +102,6 @@ class FrictionPotential : public FrictionPotentialWrapper { ...@@ -112,8 +102,6 @@ class FrictionPotential : public FrictionPotentialWrapper {
class TrivialFunction : public FrictionPotentialWrapper { class TrivialFunction : public FrictionPotentialWrapper {
public: public:
TrivialFunction() : FrictionPotentialWrapper() {}
double virtual evaluate(double) const { return 0; } double virtual evaluate(double) const { return 0; }
double virtual leftDifferential(double) const { return 0; } double virtual leftDifferential(double) const { return 0; }
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
#include "frictionpotential.hh" #include "frictionpotential.hh"
#include "localnonlinearity.hh" #include "localfriction.hh"
namespace Dune { namespace Dune {
template <class MatrixTypeTEMPLATE, class VectorTypeTEMPLATE> template <class MatrixTypeTEMPLATE, class VectorTypeTEMPLATE>
...@@ -31,7 +31,7 @@ class GlobalNonlinearity { ...@@ -31,7 +31,7 @@ class GlobalNonlinearity {
/* /*
Return a restriction of the outer function to the i'th node. Return a restriction of the outer function to the i'th node.
*/ */
virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const = 0; virtual shared_ptr<LocalFriction<dim> const> restriction(int i) const = 0;
void addHessian(VectorType const &v, MatrixType &hessian) const { void addHessian(VectorType const &v, MatrixType &hessian) const {
for (size_t i = 0; i < v.size(); ++i) { for (size_t i = 0; i < v.size(); ++i) {
......
...@@ -10,7 +10,7 @@ ...@@ -10,7 +10,7 @@
#include <dune/istl/bvector.hh> #include <dune/istl/bvector.hh>
#include "globalnonlinearity.hh" #include "globalnonlinearity.hh"
#include "localnonlinearity.hh" #include "localfriction.hh"
#include "frictionpotential.hh" #include "frictionpotential.hh"
namespace Dune { namespace Dune {
...@@ -26,12 +26,12 @@ class GlobalRuinaNonlinearity ...@@ -26,12 +26,12 @@ class GlobalRuinaNonlinearity
GlobalRuinaNonlinearity(dataref nodalIntegrals, FrictionData const &fd, GlobalRuinaNonlinearity(dataref nodalIntegrals, FrictionData const &fd,
dataref state) dataref state)
: restrictions(nodalIntegrals.size()) { : restrictions(nodalIntegrals.size()) {
auto trivialNonlinearity = make_shared<LocalNonlinearity<dim> const>( auto trivialNonlinearity = make_shared<LocalFriction<dim> const>(
make_shared<TrivialFunction const>()); make_shared<TrivialFunction const>());
for (size_t i = 0; i < restrictions.size(); ++i) { for (size_t i = 0; i < restrictions.size(); ++i) {
restrictions[i] = nodalIntegrals[i] == 0 restrictions[i] = nodalIntegrals[i] == 0
? trivialNonlinearity ? trivialNonlinearity
: make_shared<LocalNonlinearity<dim> const>( : make_shared<LocalFriction<dim> const>(
make_shared<FrictionPotential const>( make_shared<FrictionPotential const>(
nodalIntegrals[i], fd, state[i])); nodalIntegrals[i], fd, state[i]));
} }
...@@ -40,12 +40,12 @@ class GlobalRuinaNonlinearity ...@@ -40,12 +40,12 @@ class GlobalRuinaNonlinearity
/* /*
Return a restriction of the outer function to the i'th node. Return a restriction of the outer function to the i'th node.
*/ */
virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const { virtual shared_ptr<LocalFriction<dim> const> restriction(int i) const {
return restrictions[i]; return restrictions[i];
} }
private: private:
std::vector<shared_ptr<LocalNonlinearity<dim> const>> restrictions; std::vector<shared_ptr<LocalFriction<dim> const>> restrictions;
}; };
} }
#endif #endif
#ifndef DUNE_TECTONIC_LOCAL_NONLINEARITY_HH #ifndef DUNE_TECTONIC_LOCAL_FRICTION_HH
#define DUNE_TECTONIC_LOCAL_NONLINEARITY_HH #define DUNE_TECTONIC_LOCAL_FRICTION_HH
#include <cmath> #include <cmath>
#include <limits> #include <limits>
...@@ -13,15 +13,12 @@ ...@@ -13,15 +13,12 @@
#include "frictionpotential.hh" #include "frictionpotential.hh"
namespace Dune { namespace Dune {
template <int dimension> class LocalNonlinearity { template <int dimension> class LocalFriction {
public: public:
using VectorType = FieldVector<double, dimension>; using VectorType = FieldVector<double, dimension>;
using MatrixType = FieldMatrix<double, dimension, dimension>; using MatrixType = FieldMatrix<double, dimension, dimension>;
LocalNonlinearity(shared_ptr<FrictionPotentialWrapper const> func) LocalFriction(shared_ptr<FrictionPotentialWrapper const> func) : func(func) {}
: func(func) {}
std::vector<double> const &get_kinks() const { return func->get_kinks(); }
double operator()(VectorType const &x) const { double operator()(VectorType const &x) const {
return func->evaluate(x.two_norm()); return func->evaluate(x.two_norm());
......
#ifndef MINIMISATION2_HH
#define MINIMISATION2_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>
#include "minimisation.hh"
namespace Dune {
// NOTE: only works for the 2D case
template <class Functional>
void minimise2(Functional const &J, typename Functional::SmallVector &x,
size_t steps, Bisection const &bisection) {
using SmallVector = typename Functional::SmallVector;
minimisationInitialiser(J, bisection, x);
{ // Make a single step where we already know we're not differentiable
SmallVector descDir;
if (x == SmallVector(0))
J.descentAtZero(descDir);
else
J.descentDirectionNew(x, descDir);
descentMinimisation(J, x, descDir, bisection);
}
// From here on, we're in a C1 region and stay there.
for (size_t i = 1; i < steps; ++i) {
SmallVector descDir;
J.descentDirectionNew(x, descDir);
descentMinimisation(J, x, descDir, bisection);
}
}
template <class Functional>
void minimisationInitialiser(Functional const &J, Bisection const &bisection,
typename Functional::SmallVector &startingPoint) {
using SmallVector = typename Functional::SmallVector;
std::vector<double> const &kinks = J.get_kinks();
SmallVector x_old(0);
double Jx_old = J(x_old);
for (auto &kink : kinks) {
SmallVector x1 = { kink, 0 }; // Random vector that has the right norm
SmallVector const descDir1 = { x1[1], -x1[0] };
tangentialMinimisation(J, x1, descDir1, bisection);
double const Jx1 = J(x1);
SmallVector x2(0);
x2.axpy(-1, x1);
SmallVector const descDir2 = { x2[1], -x2[0] };
tangentialMinimisation(J, x2, descDir2, bisection);
double const Jx2 = J(x2);
double const Jx = (Jx2 < Jx1 ? Jx2 : Jx1);
SmallVector const x = (Jx2 < Jx1 ? x2 : x1);
if (Jx >= Jx_old)
break;
Jx_old = Jx;
x_old = x;
}
startingPoint = x_old;
}
}
#endif
...@@ -14,7 +14,6 @@ ...@@ -14,7 +14,6 @@
#include <dune/tnnmg/problem-classes/bisection.hh> #include <dune/tnnmg/problem-classes/bisection.hh>
#include "globalnonlinearity.hh" #include "globalnonlinearity.hh"
#include "localnonlinearity.hh"
#include "minimisation.hh" #include "minimisation.hh"
#include "mydirectionalconvexfunction.hh" #include "mydirectionalconvexfunction.hh"
#include "ellipticenergy.hh" #include "ellipticenergy.hh"
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment