Commit 573de679 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Move class ScalarSumFunctional to a separate file

This is a first step towards making it more general.  There is a copy
of it in ductile-fracture.cc, which I would like to eliminate.
parent 3c44697c
Pipeline #8793 passed with stage
in 5 minutes and 15 seconds
......@@ -21,6 +21,7 @@ install(FILES
integralfunctionalconstrainedlinearization.hh
miehesplitting.hh
projectednewtonstep.hh
scalarsumfunctional.hh
spectraltophysicaltransform.hh
damagedelasticenergydensity_decomposed.hh
damagedelasticenergydensity_nondecomposed.hh
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FRACTUREPHASEFIELDS_SCALARSUMFUNCTIONAL_HH
#define DUNE_FRACTUREPHASEFIELDS_SCALARSUMFUNCTIONAL_HH
#include <limits>
#include <dune/solvers/common/interval.hh>
namespace Dune::FracturePhasefields
{
// In principle this could be implemented as a simple lambda,
// but this would lead to an incomplete type warning.
template <class SSF>
class ScalarSumFunctionalDerivative
{
public:
ScalarSumFunctionalDerivative(const SSF& ssf) : ssf_(ssf) {}
double operator()(const double& z) const {
return ssf_.subDifferential(z).projectIn(0);
}
const SSF& ssf_;
};
/** \brief Sum of two functionals defined on the set of real numbers
*
* \tparam DamageFunctional Expected to model BoxConstrainedQuadraticFunctional
* \tparam DamagedElasticityFunctional Expected to be smooth but nonquadratic
*/
template <class DamageFunctional, class DamagedElasticityFunctional, class GlobalFunctional>
class ScalarSumFunctional
{
public:
ScalarSumFunctional(const DamageFunctional& damageFunctional,
const DamagedElasticityFunctional& damagedElasticityFunctional, const GlobalFunctional& J)
: damageFunctional_(damageFunctional),
damagedElasticityFunctional_(damagedElasticityFunctional),
J_(J)
{}
Solvers::Interval<double> domain() const
{
return Solvers::Interval<double>{damageFunctional_.lowerObstacle(),
damageFunctional_.upperObstacle() };
}
Solvers::Interval<double> subDifferential(double z) const
{
constexpr double tolerance = 1e-15;
// quadratic part with obstacle
Solvers::Interval<double> dJ = z*damageFunctional_.quadraticPart() - damageFunctional_.linearPart();
// Compute the subdifferential of the obstacle constraint by hand
// TODO: This is stupid -- the functional should do this for us
if (z < damageFunctional_.lowerObstacle())
dJ = -std::numeric_limits<double>::max();
if (z <= tolerance + damageFunctional_.lowerObstacle())
dJ[0] = -std::numeric_limits<double>::max();
if (z >= -tolerance + damageFunctional_.upperObstacle())
dJ[1] = std::numeric_limits<double>::max();
if (z > damageFunctional_.upperObstacle())
dJ = std::numeric_limits<double>::max();
// Smooth nonquadratic part
dJ += damagedElasticityFunctional_.subDifferential(z);
return dJ;
}
double operator()(const double& z) const
{
const auto& origin = damagedElasticityFunctional_.origin();
const auto& direction = damagedElasticityFunctional_.direction();
temp = direction;
temp *= z;
temp += origin;
return J_(temp);
}
friend auto derivative(const ScalarSumFunctional& sumFunctional)
{
// Using a lambda would lead to an incomplete type warning
return ScalarSumFunctionalDerivative(sumFunctional);
}
const DamageFunctional& damageFunctional_;
const DamagedElasticityFunctional& damagedElasticityFunctional_;
const GlobalFunctional& J_;
mutable typename DamagedElasticityFunctional::GlobalVector temp;
};
} // namespace Dune::FracturePhasefields
#endif // DUNE_FRACTUREPHASEFIELDS_SCALARSUMFUNCTIONAL_HH
......@@ -68,6 +68,7 @@
#include <dune/fracture-phasefields/damagedelasticenergydensity_decomposed.hh>
#include <dune/fracture-phasefields/damagedelasticenergydensity_nondecomposed.hh>
#include <dune/fracture-phasefields/integraldensity.hh>
#include <dune/fracture-phasefields/scalarsumfunctional.hh>
//#include <dune/matrix-vector/densematrixview.hh>
......@@ -77,99 +78,6 @@
using namespace Dune;
// In principle this could be implemented as a simple lambda,
// but this would lead to an incomplete type warning.
template <class SSF>
class ScalarSumFunctionalDerivative
{
public:
ScalarSumFunctionalDerivative(const SSF& ssf) : ssf_(ssf) {}
double operator()(const double& z) const {
return ssf_.subDifferential(z).projectIn(0);
}
const SSF& ssf_;
};
/** \brief Sum of two functionals defined on the set of real numbers
*
* \tparam DamageFunctional Expected to model BoxConstrainedQuadraticFunctional
* \tparam DamagedElasticityFunctional Expected to be smooth but nonquadratic
*/
template <class DamageFunctional, class DamagedElasticityFunctional, class GlobalFunctional>
class ScalarSumFunctional
{
public:
ScalarSumFunctional(const DamageFunctional& damageFunctional,
const DamagedElasticityFunctional& damagedElasticityFunctional, const GlobalFunctional& J)
: damageFunctional_(damageFunctional),
damagedElasticityFunctional_(damagedElasticityFunctional),
J_(J)
{}
Solvers::Interval<double> domain() const
{
return Solvers::Interval<double>{damageFunctional_.lowerObstacle(),
damageFunctional_.upperObstacle() };
}
Solvers::Interval<double> subDifferential(double z) const
{
constexpr double tolerance = 1e-15;
// quadratic part with obstacle
Solvers::Interval<double> dJ = z*damageFunctional_.quadraticPart() - damageFunctional_.linearPart();
// Compute the subdifferential of the obstacle constraint by hand
// TODO: This is stupid -- the functional should do this for us
if (z < damageFunctional_.lowerObstacle())
dJ = -std::numeric_limits<double>::max();
if (z <= tolerance + damageFunctional_.lowerObstacle())
dJ[0] = -std::numeric_limits<double>::max();
if (z >= -tolerance + damageFunctional_.upperObstacle())
dJ[1] = std::numeric_limits<double>::max();
if (z > damageFunctional_.upperObstacle())
dJ = std::numeric_limits<double>::max();
// Smooth nonquadratic part
dJ += damagedElasticityFunctional_.subDifferential(z);
return dJ;
}
double operator()(const double& z) const
{
const auto& origin = damagedElasticityFunctional_.origin();
const auto& direction = damagedElasticityFunctional_.direction();
temp = direction;
temp *= z;
temp += origin;
return J_(temp);
}
friend auto derivative(const ScalarSumFunctional& sumFunctional)
{
// Using a lambda would lead to an incomplete type warning
return ScalarSumFunctionalDerivative(sumFunctional);
}
const DamageFunctional& damageFunctional_;
const DamagedElasticityFunctional& damagedElasticityFunctional_;
const GlobalFunctional& J_;
mutable typename DamagedElasticityFunctional::GlobalVector temp;
};
class NewBisection
{
......@@ -869,7 +777,7 @@ int main (int argc, char *argv[]) try
auto bisectionSolver = [&](auto& xi, const auto& restriction, const auto& ignoreI) {
ScalarSumFunctional sumFunctional(std::get<0>(restriction.functions()), std::get<1>(restriction.functions()), J);
FracturePhasefields::ScalarSumFunctional sumFunctional(std::get<0>(restriction.functions()), std::get<1>(restriction.functions()), J);
#if 0
Bisection bisection;
......
Supports Markdown
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