Skip to content
Snippets Groups Projects
Commit 913c0611 authored by podlesny's avatar podlesny
Browse files

remove first order functional, compute eigenvalue locally in... 

remove first order functional, compute eigenvalue locally in coordinateRestriction instead of globally
parent 59547d32
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/tectonic/spatial-solving/tnnmg/functional.hh
+ 25
− 180
View file @ 913c0611
......@@ -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);
}
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment