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

Implement the NormFunctional for the new TNNMG implementation

parent 31d003bb
Pipeline #4889 failed with stage
in 9 seconds
#ifndef DUNE_TNNMG_FUNCTIONALS_NORMFUNCTIONAL_HH
#define DUNE_TNNMG_FUNCTIONALS_NORMFUNCTIONAL_HH
#include <dune/common/bitsetvector.hh>
#include <dune/solvers/common/resize.hh>
#include <dune/tnnmg/problem-classes/nonlinearity.hh>
#include <dune/tnnmg/functionals/nonsmoothconvexfunctional.hh>
namespace Dune {
namespace TNNMG {
/** \brief The weighted sum of the Euclidean norm of the vector blocks
*
* This Functional class implements the functional
* \f[
* j(v) = \sum_{i=1}^n a_i |v_i|,
* \f]
*where \f$ v \f$ is a vector in \f$ (\mathbb{R}^N)^n \f$, and \f$ a \in \mathbb{R}^n \f$
* is a vector of coefficients.
*
* \tparam blockSize The size \f$ N \f$ of the vector blocks.
*/
template <int blockSize>
class NormFunctional
: public Dune::TNNMG::NonsmoothConvexFunctional<FieldVector<double,blockSize>,
FieldMatrix<double,blockSize,blockSize> >
namespace TNNMG {
template <class V, class R=double>
class ShiftedNormFunctional;
/** \brief Coordinate restriction of the Euclidean norm functional
*
* This general implementation doesn't really allow anything but shifting.
* The shifted version be then be restricted further.
*
* \tparam V Type of the vector before restriction
* \tparam R Range type
* \tparam I Index type used to access individual vector entries
*/
template <class V, class R, class I = std::size_t>
class NormCoordinateRestriction
{
public:
using GlobalVector = V;
using GlobalIndex = I;
using Vector = std::decay_t<decltype(std::declval<GlobalVector>()[std::declval<GlobalIndex>()])>;
using Range = R;
NormCoordinateRestriction(const GlobalVector& coefficients, const GlobalVector& origin, GlobalIndex index)
: coefficients_(&coefficients),
origin_(origin),
index_(index)
{}
Range operator()(const Vector& v) const
{
DUNE_THROW(Dune::NotImplemented, "Evaluation of NormCoordinateRestriction not implemented");
}
GlobalIndex index() const
{
return index_;
}
const auto& coefficients() const
{
return (*coefficients_)[index_];
}
template<class Origin,
std::enable_if_t<! Dune::IsNumber<Origin>::value, int> = 0>
friend ShiftedNormFunctional<Vector, Range>
shift(const NormCoordinateRestriction& f, const Origin& origin)
{
return ShiftedNormFunctional<Vector, Range>(f.coefficients(),
f.origin_[f.index_], // origin of the functional before this shift
origin); // new additional shift, with reference semantics
}
protected:
const V* coefficients_;
const GlobalVector origin_;
const GlobalIndex index_;
};
/** \brief Coordinate restriction of Euclidean norm functional to a scalar function
*
* This is the specialization for the restriction to 1-dimensional subspaces.
* Unlike the general implementation above, this specialization can produce the subdifferential
* of the functional.
*
* \tparam V Type of the vector before restriction
* \tparam R Range type
* \tparam I Index type used to access individual vector entries
*
*/
template <int blocksize, class R, class I>
class NormCoordinateRestriction<FieldVector<double,blocksize>, R, I>
{
public:
using GlobalVector = FieldVector<double,blocksize>;
using GlobalIndex = I;
using Vector = std::decay_t<decltype(std::declval<GlobalVector>()[std::declval<GlobalIndex>()])>;
using Range = R;
NormCoordinateRestriction(const GlobalVector& coefficients, const GlobalVector& origin, GlobalIndex index)
: coefficients_(&coefficients),
origin_(origin),
index_(index)
{}
Range operator()(const Vector& v) const
{
auto tmp = origin_;
tmp[index_] += v;
return tmp.two_norm();
}
GlobalIndex index() const
{
return index_;
}
const auto& coefficients() const
{
return (*coefficients_)[index_];
}
void domain(Solvers::Interval<double>& d) const
{
// Our domain is the entire space
d[0] = -std::numeric_limits<double>::infinity();
d[1] = std::numeric_limits<double>::infinity();
}
void subDiff(double z, Solvers::Interval<double>& dJ) const
{
constexpr double tolerance = 1e-15;
dJ = 0.0;
// Compute norm at coefficient i
auto tmp = *origin_;
tmp[index_] += z;
double norm = tmp.two_norm();
if (norm < tolerance)
{
public:
typedef Dune::FieldVector<double,blockSize> LocalVectorType;
typedef Dune::FieldMatrix<double,blockSize,blockSize> LocalMatrixType;
typedef Nonlinearity<LocalVectorType, LocalMatrixType> NonlinearityType;
typedef typename NonlinearityType::VectorType VectorType;
typedef typename NonlinearityType::MatrixType MatrixType;
typedef typename NonlinearityType::IndexSet IndexSet;
/** \brief Constructor from a set of coefficients
*
* \param coefficients The coefficient vector \f$ a \f$
* \param tolerance The nonsmoothness a 0 of the norm function is implemented to have the radius 'tolerance'
*/
NormFunctional(const std::vector<double> coefficients,
double tolerance = 1e-15)
: coefficients_(coefficients),
tolerance_(tolerance)
{}
/** \brief Evaluate the functional at a given vector v
*/
double operator()(const VectorType& v) const
{
double result = 0.0;
// We are at 0 -- the result is a true subdifferential
dJ[0] -= (*coefficients_)[0];
dJ[1] += (*coefficients_)[0];
}
else
{
// Compute the j-th entry of the gradient
dJ += (*coefficients_)[0] * tmp[index_] / norm;
}
}
for(size_t row=0; row<v.size(); ++row)
result += coefficients_[row] * v[row].two_norm();
protected:
return result;
}
const GlobalVector* coefficients_;
/** \brief Add the gradient of the norm functional to a given vector
*
* Note: the gradient of \f$ |v| \f$ is \f$ v/|v| \f$.
*/
void addGradient(const VectorType& v, VectorType& gradient) const
{
for(size_t row=0; row<v.size(); ++row)
{
double norm = v[row].two_norm();
const GlobalVector origin_;
const GlobalIndex index_;
};
// Do nothing if we are at the nondifferentiable point
if (norm < tolerance_)
continue;
for (int j=0; j<blockSize; j++)
gradient[row][j] += coefficients_[row] * v[row][j] / norm;
}
}
/** \brief Add the Hessian of the functional to a given matrix
*
* Note: Entry \f$ i,j \f$ of the Hessian of \f$ |v| \f$ is \f[- \frac{v_i v_j}{|v|^3} + \frac{\delta_{ij}}{|v|} \f]
*/
void addHessian(const VectorType& v, MatrixType& hessian) const
{
for(size_t row=0; row<v.size(); ++row)
{
double norm = v[row].two_norm();
// Do nothing if we are at the nondifferentiable point
if (norm < tolerance_)
continue;
for (int i=0; i<blockSize; i++)
for (int j=0; j<blockSize; j++)
hessian[row][row][i][j] += coefficients_[row] * (- v[row][i] * v[row][j] / (norm*norm*norm) + (i==j) / norm);
}
}
/** \brief Shifted Euclidean norm functional
*
* \tparam V global vector type
* \tparam R Range type
*/
template <class V, class R>
class ShiftedNormFunctional
{
public:
void addHessianIndices(IndexSet& indices) const {}
using Vector = V;
using Range = R;
/** \brief Compute the subdifferential of the nonlinearity restricted to the line \f$ u + t v \f$.
*
* \param[out] subdifferential the resulting subdifferential
* \param u base point
* \param v direction
*/
virtual void directionalSubDiff(const VectorType& u, const VectorType& v, Dune::Solvers::Interval<double>& subdifferential) const
{
subdifferential[0] = 0.0;
subdifferential[1] = 0.0;
for(size_t row=0; row<u.size(); ++row)
{
double norm = u[row].two_norm();
if (norm >= tolerance_){
subdifferential[0] += coefficients_[row] * (u[row]*v[row]) / norm;
subdifferential[1] += coefficients_[row] * (u[row]*v[row]) / norm;
} else {
subdifferential[0] -= coefficients_[row] * v[row].two_norm();
subdifferential[1] += coefficients_[row] * v[row].two_norm();
}
}
}
ShiftedNormFunctional(const Vector& coefficients, const Vector& origin, const Vector& offset)
: coefficients_(&coefficients),
origin_(origin),
offset_(&offset)
{}
/** \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
*/
virtual void subDiff(int i, double x, Dune::Solvers::Interval<double>& D, int j) const
{
D[0] = 0.0;
D[1] = 0.0;
// Compute current tangential displacement
LocalVectorType displacement = u_[i];
displacement[j] = x;
double norm = displacement.two_norm();
if (norm < tolerance_)
{
D[0] = - coefficients_[i];
D[1] = coefficients_[i];
} else {
D[0] = D[1] = coefficients_[i] * displacement[j] / norm;
}
}
Range operator()(const Vector& v) const
{
DUNE_THROW(Dune::NotImplemented, "Evaluation of ShiftedNormFunctional not implemented");
}
/** \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
*/
virtual double regularity(int i, double x, int j) const
{
// Compute current tangential displacement
LocalVectorType displacement = u_[i];
displacement[j] = x;
/** \brief Tell the functional that the origin has changed */
void updateOrigin()
{
// Does not do anything, because no information is precomputed
}
return (displacement.two_norm() < tolerance_)
? std::numeric_limits<double>::max()
: 0;
}
void updateOrigin(std::size_t i)
{
// Does not do anything, because no information is precomputed
}
/** \brief Domain of definition of the functional in one coordinate direction
*
* The norm functional is the entire space
*/
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();
}
template <class Index>
friend NormCoordinateRestriction<Vector, Range, Index>
coordinateRestriction(const ShiftedNormFunctional& f, const Index& i)
{
auto tmp = f.origin_;
tmp += *f.offset_;
return NormCoordinateRestriction<Vector, Range, Index>(*f.coefficients_, tmp, i);
}
protected:
const Vector* coefficients_;
const Vector origin_;
const Vector* offset_;
};
/** \brief Directional restriction of the Euclidean norm functional
*
* \tparam V global vector type
* \tparam R Range type
*/
template<class V, class R=double>
class NormDirectionalRestriction
{
public:
using GlobalVector = V;
using Vector = typename GlobalVector::field_type;
using Range = R;
NormDirectionalRestriction(const GlobalVector& coefficients, const GlobalVector& origin, const GlobalVector& direction)
: coefficients_(&coefficients),
origin_(&origin),
direction_(&direction)
{}
Range operator()(const Vector& v) const
{
Range result(0.0);
for(size_t i=0; i<(*origin_).size(); ++i)
{
// Extract the plastic strain components into a separate data structure, for easier readability
auto p = (*origin_)[i];
p.axpy(v, (*direction_)[i]);
result += (*coefficients_)[i][0] * p.two_norm();
}
return result;
}
GlobalVector& origin()
{
return *origin_;
}
/** \brief Set the internal position vector u_pos to v.
*/
virtual void setVector(const VectorType& v)
const GlobalVector& origin() const
{
return *origin_;
}
void domain(Dune::Solvers::Interval<double>& d) const
{
d[0] = -std::numeric_limits<double>::max();
d[1] = std::numeric_limits<double>::max();
}
void subDiff(double z, Dune::Solvers::Interval<double>& dJ) const
{
double tolerance = 1e-15;
dJ = 0.0;
for (size_t row=0; row<(*origin_).size(); ++row)
{
auto p = (*origin_)[row];
p.axpy(z, (*direction_)[row]);
auto norm = p.two_norm();
if (norm >= tolerance)
{
u_ = v;
}
// Functional is differentiable here -- compute the standard directional derivative
// We compute the directional derivative as the product of the gradient and the direction
FieldVector<double,p.size()> gradient(0);
for (size_t i=0; i<p.size(); i++)
gradient[i] += (*coefficients_)[row][0] * p[i] / norm;
/** \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)
dJ += gradient * (*direction_)[row];
}
else
{
u_[i][j] = x;
dJ[0] -= (*coefficients_)[row][0] * (*direction_)[row].two_norm();
dJ[1] += (*coefficients_)[row][0] * (*direction_)[row].two_norm();
}
}
}
protected:
const GlobalVector* coefficients_;
const GlobalVector* origin_;
const GlobalVector* direction_;
};
/** \brief The weighted sum of the Euclidean norm of the vector blocks
*
* This Functional class implements the functional
* \f[
* j(v) = \sum_{i=1}^n a_i |v_i|,
* \f]
*where \f$ v \f$ is a vector in \f$ (\mathbb{R}^N)^n \f$, and \f$ a \in \mathbb{R}^n \f$
* is a vector of coefficients.
*
* \tparam V Vector type
* \tparam R Range type
*/
template<class V, class R=double>
class NormFunctional
{
public:
using Vector = V;
using Range = R;
NormFunctional(const Vector& coefficients) :
coefficients_(coefficients)
{}
Range operator()(const Vector& v) const
{
Range result(0.0);
for (size_t i=0; i<v.size(); ++i)
result += coefficients_[i][0] * v[i].two_norm();
return result;
}
private:
/** \brief For each block, the norm is weighted with the corresponding coefficient from this array */
const std::vector<double> coefficients_;
const Vector& coefficients() const
{
return coefficients_;
}
/** \brief We consider the norm function to be nondifferentiable a circle of this radius around 0 */
const double tolerance_;
friend auto directionalRestriction(const NormFunctional& f, const Vector& origin, const Vector& direction)
{
return NormDirectionalRestriction<Vector, Range>(f.coefficients_, origin, direction);
}
/** \brief The current state */
VectorType u_;
friend ShiftedNormFunctional<Vector, Range>
shift(const NormFunctional& f, const Vector& origin)
{
Vector zero;
Solvers::resizeInitialize(zero, origin, 0.0);
return ShiftedNormFunctional<Vector, Range>(f.coefficients_, std::move(zero), origin);
}
};
protected:
const Vector coefficients_;
};
}
} // namespace TNNMG
}
} // namespace Dune
#endif
......@@ -13,20 +13,26 @@ using namespace Dune;
int main(int argc, char** argv)
{
// Create a norm functional on (R^3)^2 for testing
std::vector<double> coefficients = {2.0, 3.0};
typedef TNNMG::NormFunctional<3> Functional;
Functional functional(coefficients);
const int blocksize = 3;
using Vector = BlockVector<FieldVector<double,blocksize> >;
using Weights = Vector;
// TODO: Improve: we use the same vector type for weights as for coefficient vectors,
// even though we have only one coefficient per vector block.
Weights weights = {{2.0,2.0,2.0}, {3.0,3.0,3.0}};
TNNMG::NormFunctional<Vector> functional(weights);
// A set of local test points, i.e., values for a single vector block
std::vector<Functional::LocalVectorType> localTestPoints = {{0,0,0},
{1,0,0},
{0,-2,0},
{3.14,3.14,-4}};
std::vector<FieldVector<double,blocksize> > localTestPoints = {{0,0,0},
{1,0,0},
{0,-2,0},
{3.14,3.14,-4}};
// Create real test points (i.e., block vectors) from the given values for a single block
std::vector<Functional::VectorType> testPoints(localTestPoints.size()*localTestPoints.size());
std::vector<Vector> testPoints(localTestPoints.size()*localTestPoints.size());
for (size_t i=0; i<testPoints.size(); i++)
testPoints[i].resize(coefficients.size());
testPoints[i].resize(weights.size());
for (size_t i=0; i<localTestPoints.size(); i++)
for (size_t j=0; j<localTestPoints.size(); j++)
......@@ -36,8 +42,46 @@ int main(int argc, char** argv)
}