damping.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Martin J. Kuehn, Daniel Abele
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef DAMPING_H
#define DAMPING_H
 
#include "memilio/config.h" // IWYU pragma: keep
#include "memilio/math/eigen.h" // IWYU pragma: keep
#include "memilio/utils/compiler_diagnostics.h"
#include "memilio/utils/type_safe.h"
#include "memilio/utils/stl_util.h"
#include "memilio/math/matrix_shape.h"
#include "memilio/math/smoother.h"
#include "memilio/math/floating_point.h"
 
#include <tuple>
#include <vector>
#include <algorithm>
#include <ostream>
 
namespace mio
{
 
DECL_TYPESAFE(int, DampingLevel);
 
DECL_TYPESAFE(int, DampingType);
 
template <typename FP>
class MEMILIO_ENABLE_EBO SimulationTime : public TypeSafe<FP, SimulationTime<FP>>,
                                          public OperatorAdditionSubtraction<SimulationTime<FP>>,
                                          public OperatorScalarMultiplicationDivision<SimulationTime<FP>, FP>,
                                          public OperatorComparison<SimulationTime<FP>>
{
public:
    using TypeSafe<FP, SimulationTime<FP>>::TypeSafe;
};
 
template <typename FP, class S>
class Damping : public std::tuple<typename S::Matrix, DampingLevel, DampingType, SimulationTime<FP>>
{
public:
    using Shape  = S;
    using Matrix = typename Shape::Matrix;
    using Base   = std::tuple<Matrix, DampingLevel, DampingType, SimulationTime<FP>>;
 
    template <class... T>
        requires std::is_constructible_v<Shape, T...>
    explicit Damping(T... shape_args)
        : Base(Matrix::Zero(Shape(shape_args...).rows(), Shape(shape_args...).cols()), {}, {}, {})
    {
    }
 
    template <class ME>
    Damping(const Eigen::MatrixBase<ME>& m, DampingLevel level, DampingType type, SimulationTime<FP> t)
        : Base(m, level, type, t)
    {
        assert((get_coeffs().array() <= 1.).all() && "damping coefficient out of range");
    }
 
    template <class... T>
        requires std::is_constructible_v<Shape, T...>
    Damping(FP d, DampingLevel level, DampingType type, SimulationTime<FP> t, T... shape_args)
        : Damping(Matrix::Constant(Shape(shape_args...).rows(), Shape(shape_args...).cols(), d), level, type, t)
    {
    }
 
    template <class ME>
    Damping(const Eigen::MatrixBase<ME>& m, SimulationTime<FP> t)
        : Damping(m, DampingLevel(0), DampingType(0), t)
    {
    }
 
    template <class... T>
        requires std::is_constructible_v<Shape, T...>
    Damping(FP d, SimulationTime<FP> t, T... shape_args)
        : Damping(d, DampingLevel(0), DampingType(0), t, shape_args...)
    {
    }
 
    SimulationTime<FP>& get_time()
    {
        return std::get<SimulationTime<FP>>(*this);
    }
    const SimulationTime<FP>& get_time() const
    {
        return std::get<SimulationTime<FP>>(*this);
    }
 
    DampingLevel& get_level()
    {
        return std::get<DampingLevel>(*this);
    }
    const DampingLevel& get_level() const
    {
        return std::get<DampingLevel>(*this);
    }
 
    DampingType& get_type()
    {
        return std::get<DampingType>(*this);
    }
    const DampingType& get_type() const
    {
        return std::get<DampingType>(*this);
    }
 
    const Matrix& get_coeffs() const
    {
        return std::get<Matrix>(*this);
    }
    Matrix& get_coeffs()
    {
        return std::get<Matrix>(*this);
    }
 
    Shape get_shape() const
    {
        return Shape::get_shape_of(get_coeffs());
    }
 
    friend void PrintTo(const Damping& self, std::ostream* os)
    {
        *os << '[' << std::get<SimulationTime<FP>>(self) << ',' << std::get<DampingType>(self) << ','
            << std::get<DampingLevel>(self) << ']';
        *os << '\n' << std::get<Matrix>(self);
    }
 
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("Damping");
        obj.add_element("Time", get_time());
        obj.add_element("Type", get_type());
        obj.add_element("Level", get_level());
        obj.add_element("Coeffs", get_coeffs());
    }
 
    template <class IOContext>
    static IOResult<Damping> deserialize(IOContext& io)
    {
        auto obj = io.expect_object("Damping");
        auto ti  = obj.expect_element("Time", Tag<SimulationTime<FP>>{});
        auto ty  = obj.expect_element("Type", Tag<DampingType>{});
        auto l   = obj.expect_element("Level", Tag<DampingLevel>{});
        auto c   = obj.expect_element("Coeffs", Tag<Matrix>{});
        return apply(
            io,
            [](auto&& ti_, auto&& ty_, auto&& l_, auto&& c_) {
                return Damping(c_, l_, ty_, ti_);
            },
            ti, ty, l, c);
    }
};
 
template <typename FP, class D>
class Dampings
{
public:
    using Shape           = typename D::Shape;
    using Matrix          = typename Shape::Matrix;
    using value_type      = D;
    using reference       = value_type&;
    using const_reference = const value_type&;
    using iterator        = typename std::vector<value_type>::iterator;
    using const_iterator  = typename std::vector<value_type>::const_iterator;
 
    template <class... T>
        requires std::is_constructible_v<Shape, T...>
    explicit Dampings(T... shape_args)
        : m_dampings()
        , m_shape(shape_args...)
    {
        update_cache();
    }
 
    Dampings(std::initializer_list<value_type> il)
        : m_dampings(il)
    {
        assert(il.size() > 0);
        m_shape = il.begin()->get_shape();
        update_cache();
    }
 
    void add(const value_type& damping)
    {
        add_(damping);
    }
    template <class... T>
    void add(T&&... t)
    {
        add_(value_type(std::forward<T>(t)...));
    }
    template <class... T>
    void add(FP d, T... t)
    {
        add_(value_type(d, std::forward<T>(t)..., m_shape));
    }
 
    void remove(size_t i)
    {
        assert(m_dampings.size() > i);
        m_dampings.erase(m_dampings.begin() + i);
        automatic_cache_update();
    }
 
    void clear()
    {
        m_dampings.clear();
        automatic_cache_update();
    }
 
    void set_automatic_cache_update(bool b)
    {
        m_automatic_cache_update = b;
        automatic_cache_update();
    }
 
    auto get_matrix_at(SimulationTime<FP> t) const
    {
        assert(!m_accumulated_dampings_cached.empty() &&
               "Cache is not current. Did you disable the automatic cache update?");
 
        auto ub = std::upper_bound(m_accumulated_dampings_cached.begin(), m_accumulated_dampings_cached.end(),
                                   std::make_tuple(t), [](auto&& tup1, auto&& tup2) {
                                       return std::get<SimulationTime<FP>>(tup1) < std::get<SimulationTime<FP>>(tup2);
                                   });
 
        auto damping = smoother_cosine<FP>(
            t.get(), (std::get<SimulationTime<FP>>(*ub) - mio::SimulationTime<FP>(1.0)).get(),
            std::get<SimulationTime<FP>>(*ub).get(), std::get<Matrix>(*(ub - 1)), std::get<Matrix>(*ub));
 
        return damping;
    }
 
    reference operator[](size_t i)
    {
        return m_dampings[i];
    }
    const_reference operator[](size_t i) const
    {
        return m_dampings[i];
    }
 
    bool operator==(const Dampings& other) const
    {
        return m_dampings == other.m_dampings;
    }
    bool operator!=(const Dampings& other) const
    {
        return !(*this == other);
    }
 
    size_t get_num_dampings() const
    {
        return m_dampings.size();
    }
 
    Shape get_shape() const
    {
        return m_shape;
    }
 
    iterator begin()
    {
        return m_dampings.begin();
    }
    iterator end()
    {
        return m_dampings.end();
    }
    const_iterator begin() const
    {
        return m_dampings.begin();
    }
    const_iterator end() const
    {
        return m_dampings.end();
    }
 
    friend void PrintTo(const Dampings& self, std::ostream* os)
    {
        for (auto& d : self.m_dampings) {
            *os << '\n'
                << '[' << std::get<SimulationTime<FP>>(d) << ',' << std::get<DampingType>(d) << ','
                << std::get<DampingLevel>(d) << ']';
            *os << '\n' << std::get<Matrix>(d);
        }
    }
 
    template <class IOContext>
    void serialize(IOContext& io) const
    {
        auto obj = io.create_object("Dampings");
        obj.add_element("Shape", get_shape());
        obj.add_list("Dampings", begin(), end());
    }
 
    template <class IOContext>
    static IOResult<Dampings> deserialize(IOContext& io)
    {
        auto obj = io.expect_object("Dampings");
        auto s   = obj.expect_element("Shape", Tag<Shape>{});
        auto d   = obj.expect_list("Dampings", Tag<value_type>{});
        return apply(
            io,
            [](auto&& s_, auto&& d_) -> IOResult<Dampings> {
                Dampings dampings(s_);
                for (auto&& i : d_) {
                    if (i.get_shape() != s_) {
                        return failure(StatusCode::InvalidValue, "Dampings must all have the same shape.");
                    }
                    dampings.add(i);
                }
                return success(dampings);
            },
            s, d);
    }
 
private:
    void add_(const value_type& damping);
 
    void update_cache();
 
    void automatic_cache_update()
    {
        if (m_automatic_cache_update) {
            update_cache();
        }
    }
 
    static void update_active_dampings(
        const value_type& damping,
        std::vector<std::tuple<std::reference_wrapper<const Matrix>, DampingLevel, DampingType>>& active_by_type,
        std::vector<std::tuple<Matrix, DampingLevel>>& sum_by_level);
 
    template <class Iter>
    static void inclusive_exclusive_sum_rec(Iter b, Iter e, Matrix& sum)
    {
        if (b != e) {
            auto& mat_b   = std::get<Matrix>(*b);
            auto mat_prod = (sum.array() * mat_b.array()).matrix();
            sum           = sum + mat_b - mat_prod;
            inclusive_exclusive_sum_rec(++b, e, sum);
        }
    }
    template <class Tuple>
    static Matrix inclusive_exclusive_sum(const std::vector<Tuple>& v)
    {
        assert(!v.empty());
        Matrix sum = std::get<Matrix>(v.front());
        inclusive_exclusive_sum_rec(v.begin() + 1, v.end(), sum);
        return sum;
    }
 
    std::vector<value_type> m_dampings;
    Shape m_shape;
    std::vector<std::tuple<Matrix, SimulationTime<FP>>> m_accumulated_dampings_cached;
    bool m_automatic_cache_update = true;
};
 
template <typename FP, class D>
void Dampings<FP, D>::update_cache()
{
    using std::get;
 
    if (m_accumulated_dampings_cached.empty()) {
        m_accumulated_dampings_cached.emplace_back(Matrix::Zero(m_shape.rows(), m_shape.cols()),
                                                   SimulationTime<FP>(std::numeric_limits<FP>::lowest()));
 
        std::vector<std::tuple<std::reference_wrapper<const Matrix>, DampingLevel, DampingType>> active_by_type;
        std::vector<std::tuple<Matrix, DampingLevel>> sum_by_level;
        for (auto& damping : m_dampings) {
            //update active damping
            update_active_dampings(damping, active_by_type, sum_by_level);
            auto combined_damping = inclusive_exclusive_sum(sum_by_level);
            assert((combined_damping.array() <= FP(1.0)).all() &&
                   "unexpected error, accumulated damping out of range.");
            if (floating_point_equal<FP>(get<SimulationTime<FP>>(damping).get(),
                                         get<SimulationTime<FP>>(m_accumulated_dampings_cached.back()).get(), 1e-15,
                                         1e-15)) {
                get<Matrix>(m_accumulated_dampings_cached.back()) = combined_damping;
            }
            else {
                m_accumulated_dampings_cached.emplace_back(combined_damping, get<SimulationTime<FP>>(damping));
            }
        }
 
        m_accumulated_dampings_cached.emplace_back(get<Matrix>(m_accumulated_dampings_cached.back()),
                                                   SimulationTime<FP>(std::numeric_limits<FP>::max()));
    }
}
 
template <typename FP, class D>
void Dampings<FP, D>::add_(const value_type& damping)
{
    assert(damping.get_shape() == m_shape && "Inconsistent matrix shape.");
    insert_sorted_replace<value_type>(m_dampings, damping, [](auto& tup1, auto& tup2) {
        return std::make_tuple(tup1.get_time(), int(tup1.get_type()), int(tup1.get_level())) <
               std::make_tuple(tup2.get_time(), int(tup2.get_type()), int(tup2.get_level()));
    });
    m_accumulated_dampings_cached.clear();
    if (m_automatic_cache_update) {
        update_cache();
    }
}
 
template <typename FP, class S>
void Dampings<FP, S>::update_active_dampings(
    const value_type& damping,
    std::vector<std::tuple<std::reference_wrapper<const Matrix>, DampingLevel, DampingType>>& active_by_type,
    std::vector<std::tuple<Matrix, DampingLevel>>& sum_by_level)
{
    using std::get;
 
    const int MatrixIdx = 0;
 
    //find active with same type and level if existent
    auto iter_active_same_type = std::find_if(active_by_type.begin(), active_by_type.end(), [&damping](auto& active) {
        return get<DampingLevel>(active) == get<DampingLevel>(damping) &&
               get<DampingType>(active) == get<DampingType>(damping);
    });
    if (iter_active_same_type != active_by_type.end()) {
        //replace active of the same type and level
        auto& active_same_type = *iter_active_same_type;
        //find active with the same level
        auto& sum_same_level = *std::find_if(sum_by_level.begin(), sum_by_level.end(), [&damping](auto& sum) {
            return get<DampingLevel>(sum) == get<DampingLevel>(damping);
        });
        //remove active with the same type and level and add new one
        get<MatrixIdx>(sum_same_level) += get<MatrixIdx>(damping) - get<MatrixIdx>(active_same_type).get();
        get<MatrixIdx>(active_same_type) = get<MatrixIdx>(damping);
    }
    else {
        //add new type to active
        active_by_type.emplace_back(get<MatrixIdx>(damping), get<DampingLevel>(damping), get<DampingType>(damping));
        //find damping with same level if existent
        auto iter_sum_same_level = std::find_if(sum_by_level.begin(), sum_by_level.end(), [&damping](auto& sum) {
            return get<DampingLevel>(sum) == get<DampingLevel>(damping);
        });
        if (iter_sum_same_level != sum_by_level.end()) {
            //add to existing level
            get<MatrixIdx>(*iter_sum_same_level) += get<MatrixIdx>(damping);
        }
        else {
            //add new level
            sum_by_level.emplace_back(get<MatrixIdx>(damping), get<DampingLevel>(damping));
        }
    }
}
 
template <typename FP>
using SquareDamping = Damping<FP, SquareMatrixShape<FP>>;
template <typename FP>
using SquareDampings = Dampings<FP, SquareDamping<FP>>;
template <typename FP>
using VectorDamping = Damping<FP, ColumnVectorShape<FP>>;
template <typename FP>
using VectorDampings = Dampings<FP, VectorDamping<FP>>;
 
} // namespace mio
 
#endif // DAMPING_H