Skip to content

GitLab

Commit d839b937 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Remove the TrescaFriction implementation for the old interface 

It is very old, after all.  Unfortunately, it is the one that
is described in the documentation, but that documentation
will need to be rewritten from scratch anyway.
parent 9a611b3e
Hide whitespace changes
Inline Side-by-side
dune/tnnmg/functionals/trescafrictionfunctional.hh
View file @ d839b937
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=8 sw=2 et sts=2:
#ifndef DUNE_TNNMG_EXAMPLES_TRESCAFRICTIONFUNCTIONAL_HH
#define DUNE_TNNMG_EXAMPLES_TRESCAFRICTIONFUNCTIONAL_HH
#ifndef DUNE_TNNMG_FUNCTIONALS_TRESCAFRICTIONFUNCTIONAL_HH
#define DUNE_TNNMG_FUNCTIONALS_TRESCAFRICTIONFUNCTIONAL_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/typeutilities.hh>
......@@ -14,11 +14,6 @@
#include <dune/tnnmg/functionals/nonsmoothconvexfunctional.hh>
#ifndef NEW_TRESCA_FRICTION_FUNCTIONAL
#include <dune/tnnmg/problem-classes/nonlinearity.hh>
#endif
#ifdef NEW_TRESCA_FRICTION_FUNCTIONAL
template <class V, class R=double>
class ShiftedTrescaFrictionFunctional;
......@@ -379,261 +374,6 @@ public:
protected:
const Vector coefficients_;
};
#else
// { trescafunctional_header_begin }
template <int dim>
class TrescaFrictionFunctional
: public Dune::TNNMG::NonsmoothConvexFunctional<Dune::FieldVector<double,dim>, Dune::FieldMatrix<double,dim,dim> >
// { trescafunctional_header_end }
{
public:
typedef Dune::FieldVector<double,dim> LocalVectorType;
typedef Dune::FieldMatrix<double,dim,dim> LocalMatrixType;
private:
// { trescafunctional_tangentialpart_begin }
static Dune::FieldVector<double,dim-1> tangentialPart(const LocalVectorType& v)
{
Dune::FieldVector<double,dim-1> result;
std::copy(v.begin()+1,v.end(),result.begin());
return result;
}
// { trescafunctional_tangentialpart_end }
public:
typedef Nonlinearity<LocalVectorType, LocalMatrixType> NonlinearityType;
typedef typename NonlinearityType::VectorType VectorType;
typedef typename NonlinearityType::MatrixType MatrixType;
typedef typename NonlinearityType::IndexSet IndexSet;
TrescaFrictionFunctional(const std::vector<double> coefficients,
const Dune::BitSetVector<1>& frictionNodes,
double tolerance = 1e-15) :
coefficients_(coefficients),
tolerance_(tolerance),
frictionNodes_(frictionNodes)
{}
// { trescafunctional_energy_begin }
double operator()(const VectorType& v) const
{
double result = 0.0;
for(size_t row=0; row<v.size(); ++row)
if (frictionNodes_[row][0])
result += coefficients_[row] * tangentialPart(v[row]).two_norm();
return result;
}
// { trescafunctional_energy_end }
/** \brief Add the gradient of the Tresca friction potential to a given vector
*
* Note: the gradient of |v| is v/|v|
*/
// { trescafunctional_addgradient_begin }
void addGradient(const VectorType& v, VectorType& gradient) const
{
for(size_t row=0; row<v.size(); ++row)
{
if (not frictionNodes_[row][0])
continue;
double tangentialNorm = tangentialPart(v[row]).two_norm();
// Do nothing if we are at a nondifferentiable point
if (tangentialNorm < tolerance_)
continue;
for (int j=1; j<dim; j++)
gradient[row][j] += coefficients_[row] * v[row][j] / tangentialNorm;
}
}
// { trescafunctional_addgradient_end }
/** \brief Add the Hessian of the Tresca friction potential to a given matrix
*
* Note: Entry ij of the Hessian of |v| is - \frac{v_i v_j}{\norm{v}^3} + \frac{\delta_{ij}}{\norm{v}}
*/
// { trescafunctional_addhessian_begin }
void addHessian(const VectorType& v, MatrixType& hessian) const
{
for(size_t row=0; row<v.size(); ++row)
{
if (not frictionNodes_[row][0])
continue;
double tangentialNorm = tangentialPart(v[row]).two_norm();
// Do nothing if we are at the nondifferentiable point
if (tangentialNorm < tolerance_)
continue;
for (int i=1; i<dim; i++)
for (int j=1; j<dim; j++)
hessian[row][row][i][j] += coefficients_[row] * (- v[row][i] * v[row][j] / (tangentialNorm*tangentialNorm*tangentialNorm) + (i==j) / tangentialNorm);
}
}
// { trescafunctional_addhessian_end }
void addHessianIndices(IndexSet& indices) const {}
/** \brief Compute the subdifferential of the nonlinearity restricted to the line u + t v.
*
* \param[out] subdifferential the resulting subdifferential
* \param u base point
* \param v direction
*/
// { trescafunctional_directional_subdiff_begin }
virtual void directionalSubDiff(const VectorType& u, const VectorType& v, Dune::Solvers::Interval<double>& subdifferential)
{
subdifferential[0] = 0.0;
subdifferential[1] = 0.0;
for(size_t row=0; row<u.size(); ++row)
{
if (not frictionNodes_[row][0])
continue;
double tangentialNorm = tangentialPart(u[row]).two_norm();
if (tangentialNorm >= tolerance_){
for (int j=1; j<dim; j++) {
subdifferential[0] += coefficients_[row] * (u[row][j]*v[row][j]) / tangentialNorm;
subdifferential[1] += coefficients_[row] * (u[row][j]*v[row][j]) / tangentialNorm;
}
} else {
for (int j=1; j<dim; j++) {
subdifferential[0] -= coefficients_[row] * v[row][j];
subdifferential[1] += coefficients_[row] * v[row][j];
}
}
}
}
// { trescafunctional_directional_subdiff_end }
/** \brief Compute the subdifferential of the nonlinearity restricted to the
* line u_pos' +t e_(i,j) at t=0.
*
* Here e_(i,j) is the (i,j)-th Euclidean unit vector,
* and u_pos' is the internal position vector u_pos with the (i,j)-the entry replaced by x.
* If the nonlinearity decouples in the Euclidean directions this is simply the (i,j)-th
* component of the subdifferential.
*
* \param i global index
* \param j local index
* \param x value of the (i,j)-th entry of position to evaluate the nonlinearity at
* \param[out] D the subdifferential
*/
// { trescafunctional_subdiff_begin }
virtual void subDiff(int i, double x, Dune::Solvers::Interval<double>& D, int j) const
{
D[0] = 0.0;
D[1] = 0.0;
if (not frictionNodes_[i][0])
return;
if (j!=0)
{
// Compute current tangential displacement
LocalVectorType displacement = u_[i];
displacement[j] = x;
double tangentialNorm = tangentialPart(displacement).two_norm();
if (tangentialNorm < tolerance_)
{
D[0] = - coefficients_[i];
D[1] = coefficients_[i];
} else
D[0] = D[1] = coefficients_[i] * displacement[j] / tangentialNorm;
}
}
// { trescafunctional_subdiff_end }
/** \brief Return the regularity of the nonlinearity restricted to the
* line u_pos' +t e_(i,j) at t=0.
*
* Here e_(i,j) is the (i,j)-th Euclidean unit vector,
* and u_pos' is the internal position vector u_pos with the (i,j)-the entry replaced by x.
* Usually this will be the third derivative or a local Lipschitz constant of the second
* derivative. Note that if the subdifferential is set-valued at this position, this
* value will normally be infinity.
*
* \param i global index
* \param j local index
* \param x value of the (i,j)-th entry of position to evaluate the nonlinearity at
* \returns a value measuring the regularity
*/
// { trescafunctional_regularity_begin }
virtual double regularity(int i, double x, int j) const
{
if (not frictionNodes_[i][0])
return 0.0;
// Compute current tangential displacement
LocalVectorType displacement = u_[i];
displacement[j] = x;
double tangentialNorm = tangentialPart(displacement).two_norm();
if (tangentialNorm < tolerance_)
return std::numeric_limits<double>::max();
return 0.0;
}
// { trescafunctional_regularity_end }
// { trescafunctional_domain_begin }
void domain(int i, Dune::Solvers::Interval<double>& dom, int j) const
{
// Our domain is the entire space
dom[0] = -std::numeric_limits<double>::max();
dom[1] = std::numeric_limits<double>::max();
}
// { trescafunctional_domain_end }
/** \brief Set the internal position vector u_pos to v.
*
* This method is empty, because it is only needed if the nonlinearity
* does not decouple in the Euclidean directions.
*/
// { trescafunctional_state_begin }
virtual void setVector(const VectorType& v)
{
u_ = v;
};
/** \brief Update the (i,j)-th entry of the internal position vector u_pos to x.
*
*
* \param i global index
* \param j local index
* \param x new value of the entry (u_pos)_(i,j)
*/
virtual void updateEntry(int i, double x, int j)
{
u_[i][j] = x;
};
// { trescafunctional_state_end }
private:
/** \brief For each block, the norm is weighted with the corresponding coefficient from this array */
const std::vector<double> coefficients_;
/** \brief We consider the norm function to be nondifferentiable a circle of this radius around 0 */
const double tolerance_;
/** \brief Is 'true' for all degrees of freedom with friction */
Dune::BitSetVector<1> frictionNodes_;
/** \brief The current state */
VectorType u_;
};
#endif
#endif
dune/tnnmg/test/trescafrictionfunctionaltest.cc
View file @ d839b937
#include "config.h"
#define NEW_TRESCA_FRICTION_FUNCTIONAL
#include <iostream>
#include <dune/common/exceptions.hh>
......
examples/trescafriction.cc
View file @ d839b937
......@@ -2,8 +2,6 @@
#include <array>
#define NEW_TRESCA_FRICTION_FUNCTIONAL
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/hybridutilities.hh>
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment