diff --git a/dune/tectonic/spatial-solving/tnnmg/functional.hh b/dune/tectonic/spatial-solving/tnnmg/functional.hh
index 15674fb378df4337709f4f3244f04a69cf76765e..4613c3aa533e8fa20ec730af569bb11475e2d5f0 100755
--- a/dune/tectonic/spatial-solving/tnnmg/functional.hh
+++ b/dune/tectonic/spatial-solving/tnnmg/functional.hh
@@ -22,7 +22,7 @@
template<class M, class V, class N, class L, class U, class R>
class DirectionalRestriction;
-template<class M, class E, class V, class N, class L, class U, class O, class R>
+template<class M, class V, class N, class L, class U, class O, class R>
class ShiftedFunctional;
template <class M, class V, class N, class L, class U, class R>
@@ -36,12 +36,11 @@ class Functional;
* \tparam V nonlinearity type
* \tparam R Range type
*/
-template<class M, class E, class V, class N, class L, class U, class O, class R>
+template<class M, class V, class N, class L, class U, class O, class R>
class ShiftedFunctional {
public:
using Matrix = std::decay_t<M>;
- using Eigenvalues = E;
using Vector = std::decay_t<V>;
using Nonlinearity = std::decay_t<N>;
using LowerObstacle = std::decay_t<L>;
@@ -50,10 +49,9 @@ class ShiftedFunctional {
using Range = R;
template <class MM, class VV, class LL, class UU, class OO>
- ShiftedFunctional(MM&& matrix, const Eigenvalues& maxEigenvalues, VV&& linearPart, const Nonlinearity& phi,
+ ShiftedFunctional(MM&& matrix, VV&& linearPart, const Nonlinearity& phi,
LL&& lower, UU&& upper, OO&& origin) :
quadraticPart_(std::forward<MM>(matrix)),
- maxEigenvalues_(maxEigenvalues),
originalLinearPart_(std::forward<VV>(linearPart)),
phi_(phi),
originalLowerObstacle_(std::forward<LL>(lower)),
@@ -118,13 +116,8 @@ class ShiftedFunctional {
return phi_;
}
- const auto& maxEigenvalues() const {
- return maxEigenvalues_;
- }
-
protected:
Dune::Solvers::ConstCopyOrReference<M> quadraticPart_;
- const Eigenvalues& maxEigenvalues_;
Dune::Solvers::ConstCopyOrReference<V> originalLinearPart_;
const Nonlinearity& phi_;
@@ -222,126 +215,6 @@ class DirectionalRestriction :
};
-/** \brief Box constrained quadratic functional with nonlinearity
- * Note: call setIgnore() to set up functional in affine subspace with ignore information
- *
- * \tparam V Vector type
- * \tparam N Nonlinearity type
- * \tparam L Lower obstacle type
- * \tparam U Upper obstacle type
- * \tparam R Range type
- */
-template<class N, class V, class R>
-class FirstOrderFunctional {
- using Interval = typename Dune::Solvers::Interval<R>;
-
-public:
- using Nonlinearity = std::decay_t<N>;
- using Vector = V;
- using LowerObstacle = V;
- using UpperObstacle = V;
- using Range = R;
-
- using field_type = typename V::field_type;
-public:
- template <class LL, class UU, class OO, class DD>
- FirstOrderFunctional(
- const Range& maxEigenvalue,
- const Range& linearPart,
- const Nonlinearity& phi,
- LL&& lower,
- UU&& upper,
- OO&& origin,
- DD&& direction) :
- quadraticPart_(maxEigenvalue),
- linearPart_(linearPart),
- lower_(std::forward<LL>(lower)),
- upper_(std::forward<UU>(upper)),
- origin_(std::forward<OO>(origin)),
- direction_(std::forward<DD>(direction)),
- phi_(phi) {
-
- // set defect obstacles
- defectLower_ = -std::numeric_limits<field_type>::max();
- defectUpper_ = std::numeric_limits<field_type>::max();
- Dune::TNNMG::directionalObstacles(origin_.get(), direction_.get(), lower_.get(), upper_.get(), defectLower_, defectUpper_);
-
- // set domain
- phi_.directionalDomain(origin_.get(), direction_.get(), domain_);
-
- if (domain_[0] < defectLower_) {
- domain_[0] = defectLower_;
- }
-
- if (domain_[1] > defectUpper_) {
- domain_[1] = defectUpper_;
- }
- }
-
- Range operator()(const Vector& v) const
- {
- DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderFunctional not implemented");
- }
-
- auto subDifferential(double x) const {
- Interval Di;
- Vector uxv = origin_.get();
-
- Dune::MatrixVector::addProduct(uxv, x, direction_.get());
- phi_.directionalSubDiff(uxv, direction_.get(), Di);
-
- const auto Axmb = quadraticPart_ * x - linearPart_;
- Di[0] += Axmb;
- Di[1] += Axmb;
-
- return Di;
- }
-
- const Interval& domain() const {
- return domain_;
- }
-
- const auto& origin() const {
- return origin_.get();
- }
-
- const auto& direction() const {
- return direction_.get();
- }
-
- const auto& lowerObstacle() const {
- return defectLower_;
- }
-
- const auto& upperObstacle() const {
- return defectUpper_;
- }
-
- const auto& quadraticPart() const {
- return quadraticPart_;
- }
-
- const auto& linearPart() const {
- return linearPart_;
- }
-private:
- const Range quadraticPart_;
- const Range linearPart_;
-
- Dune::Solvers::ConstCopyOrReference<LowerObstacle> lower_;
- Dune::Solvers::ConstCopyOrReference<UpperObstacle> upper_;
-
- Dune::Solvers::ConstCopyOrReference<Vector> origin_;
- Dune::Solvers::ConstCopyOrReference<Vector> direction_;
-
- const Nonlinearity& phi_;
-
- Interval domain_;
-
- 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.
@@ -358,20 +231,26 @@ auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
using namespace Dune::MatrixVector;
namespace H = Dune::Hybrid;
- const LocalMatrix* Aii_p = nullptr;
-
//print(f.originalLinearPart(), "f.linearPart: ");
//print(f.quadraticPart(), "f.quadraticPart: ");
+ const auto& origini = f.origin()[i];
LocalVector ri = f.originalLinearPart()[i];
+
const auto& Ai = f.quadraticPart()[i];
+
+ auto Aii = Ai[i];
+ auto maxEig = maxEigenvalue(Aii);
+
sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
- // TODO Here we must implement a wrapper to guarantee that this will work with proxy matrices!
- H::ifElse(H::equals(j, i), [&](auto&& id){
- Aii_p = id(&Aij);
- });
Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
});
+ Dune::MatrixVector::addProduct(ri, maxEig, origini);
+
+ Aii = 0;
+ for (size_t j=0; j<Aii.N(); j++) {
+ Aii[j] = maxEig;
+ }
//print(*Aii_p, "Aii_p:");
//print(ri, "ri:");
@@ -381,7 +260,12 @@ auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
auto& phii = f.phi().restriction(i);
- return ShiftedFunctional<LocalMatrix&, Range, LocalVector, std::decay_t<decltype(phii)>, LocalLowerObstacle&, LocalUpperObstacle&, LocalVector&, Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), phii, f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
+ LocalLowerObstacle dli = f.originalLowerObstacle()[i];
+ LocalUpperObstacle dui = f.originalUpperObstacle()[i];
+ dli -= origini;
+ dui -= origini;
+
+ return Functional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui));
}
@@ -399,7 +283,6 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
private:
using Base = Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>;
-
public:
using Nonlinearity = std::decay_t<N>;
@@ -409,11 +292,8 @@ 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:
@@ -425,47 +305,13 @@ 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)),
- maxEigenvalues_(this->linearPart().size())
- {
- for (size_t i=0; i<this->quadraticPart().N(); ++i) {
- const auto& Aii = this->quadraticPart()[i][i];
- maxEigenvalues_[i] = maxEigenvalue(Aii);
-
- // due to numerical imprecision, Aii might not be symmetric leading to complex eigenvalues
- // consider Toeplitz decomposition of Aii and take only symmetric part
-
- /*auto symBlock = Aii;
-
- for (size_t j=0; j<symBlock.N(); j++) {
- for (size_t k=0; k<symBlock.M(); k++) {
- if (symBlock[j][k]/symBlock[j][j] < 1e-14)
- symBlock[j][k] = 0;
- }
- }
-
- try {
- typename Vector::block_type eigenvalues;
- Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
- maxEigenvalues_[i] =
- *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
- } catch (Dune::Exception &e) {
- print(symBlock, "symBlock");
- typename Vector::block_type eigenvalues;
- Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
- std::cout << "complex eig in dof: " << i << std::endl;
- } */
- }
- }
+ phi_(std::forward<NN>(phi))
+ {}
const Nonlinearity& phi() const {
return phi_.get();
}
- const auto& maxEigenvalues() const {
- return maxEigenvalues_;
- }
-
Range operator()(const Vector& v) const
{
//print(v, "v:");
@@ -490,14 +336,13 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
dir = 0.0;
}
-
return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, dir, dirNorm);
}
friend auto shift(const Functional& f, const Vector& origin)
- -> ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>
+ -> ShiftedFunctional<Matrix&, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>
{
- return ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+ return ShiftedFunctional<Matrix&, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
}
};