b2newmark.H

//------------------------------------------------------------------------
// b2newmark.H --
//
//
// written by Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//          Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code. The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef B2NEWMARK_H_
#define B2NEWMARK_H_
 
#include "b2time_integrator.H"
 
namespace b2000::solver {
 
template <typename T>
class NewmarkIntegrator : public TimeIntegrator<T> {
public:
    using type_t = ObjectTypeComplete<NewmarkIntegrator, typename TimeIntegrator<T>::type_t>;
 
    virtual void SetAttribute(const Case& case__) override{/*TODO*/};
 
    virtual void InitializeTimeIntegrator(
          const b2linalg::Matrix<T, b2linalg::Mrectangle>& initial_dof, const double dt) override;
 
    virtual void InitializeFieldsForThisTimeStep(
          b2linalg::Matrix<T, b2linalg::Mrectangle>& dof,
          const b2linalg::Matrix<T, b2linalg::Mrectangle>& dU) override;
 
    virtual void UpdateFields(
          b2linalg::Matrix<T, b2linalg::Mrectangle>& dof,
          const b2linalg::Matrix<T, b2linalg::Mrectangle>& dU) override;
 
    virtual TimeIntegrator<T>::EffectiveSymmetricSystem GetEffectiveElementSystem(
          const b2linalg::Matrix<T, b2linalg::Mpacked>& k,
          const b2linalg::Matrix<T, b2linalg::Mpacked>& d,
          const b2linalg::Matrix<T, b2linalg::Mpacked>& m,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
          const b2linalg::Index& index) const override;
 
    virtual TimeIntegrator<T>::EffectiveUnsymmetricSystem GetEffectiveElementSystem(
          const b2linalg::Matrix<T, b2linalg::Mrectangle>& k,
          const b2linalg::Matrix<T, b2linalg::Mrectangle>& d,
          const b2linalg::Matrix<T, b2linalg::Mrectangle>& m,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
          const b2linalg::Index& index) const override;
 
    static type_t type;
 
private:
    T beta{0.25};  
    T gamma{0.5};  
};
 
template <typename T>
inline void NewmarkIntegrator<T>::InitializeTimeIntegrator(
      const b2linalg::Matrix<T, b2linalg::Mrectangle>& initial_dof, const double dt) {
    this->dt = dt;
 
    this->V_.resize(initial_dof.size1(), 1);
    this->Z_.resize(initial_dof.size1(), 1);
    this->dot_Z_.resize(initial_dof.size1(), 1);
}
 
template <typename T>
inline void NewmarkIntegrator<T>::InitializeFieldsForThisTimeStep(
      b2linalg::Matrix<T, b2linalg::Mrectangle>& dof,
      const b2linalg::Matrix<T, b2linalg::Mrectangle>& dU) {
    this->V_[0] = dof[0];      
    this->Z_[0] = dof[1];      
    this->dot_Z_[0] = dof[2];  
 
    UpdateFields(dof, dU);
}
 
template <typename T>
inline void NewmarkIntegrator<T>::UpdateFields(
      b2linalg::Matrix<T, b2linalg::Mrectangle>& dof,
      const b2linalg::Matrix<T, b2linalg::Mrectangle>& dU) {
    dof[1] = gamma / beta / this->dt * dU[0] + (1.0 - gamma / beta) * this->Z_[0];
    dof[1] += (1.0 - gamma / 2.0 / beta) * this->dt * this->dot_Z_[0];
 
    dof[2] = 1.0 / beta / this->dt / this->dt * dU[0]
             - (1.0 - 2 * beta) / 2.0 / beta * this->dot_Z_[0];
    dof[2] -= 1.0 / beta / this->dt * this->Z_[0];
}
 
template <typename T>
inline TimeIntegrator<T>::EffectiveSymmetricSystem NewmarkIntegrator<T>::GetEffectiveElementSystem(
      const b2linalg::Matrix<T, b2linalg::Mpacked>& k,
      const b2linalg::Matrix<T, b2linalg::Mpacked>& d,
      const b2linalg::Matrix<T, b2linalg::Mpacked>& m,
      const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
      const b2linalg::Index& index) const {
    b2linalg::Matrix<T, b2linalg::Mpacked> k_eff{m.size1(), m.is_upper()};
    b2linalg::Vector<T, b2linalg::Vdense> f_eff_dyn{m.size1()};
 
    k_eff = 1.0 / beta / this->dt / this->dt * m;                  
    if (d.size1() != 0) { k_eff += gamma / beta / this->dt * d; }  
    k_eff += k;                                                    
 
    b2linalg::Vector<T, b2linalg::Vdense> a{index.size()};
    for (size_t i{}; i < index.size(); ++i) { a[i] = dof(index[i], 2); }
 
    f_eff_dyn = m * a;
 
    return std::make_pair(k_eff, f_eff_dyn);
}
 
template <typename T>
inline TimeIntegrator<T>::EffectiveUnsymmetricSystem
      NewmarkIntegrator<T>::GetEffectiveElementSystem(
            const b2linalg::Matrix<T, b2linalg::Mrectangle>& k,
            const b2linalg::Matrix<T, b2linalg::Mrectangle>& d,
            const b2linalg::Matrix<T, b2linalg::Mrectangle>& m,
            const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
            const b2linalg::Index& index) const {
    b2linalg::Matrix<T, b2linalg::Mrectangle> k_eff{m.size1(), m.size2()};
    b2linalg::Vector<T, b2linalg::Vdense> f_eff_dyn{m.size1()};
 
    k_eff = 1.0 / beta / this->dt / this->dt * m;                  
    if (d.size1() != 0) { k_eff += gamma / beta / this->dt * d; }  
    k_eff += k;                                                    
 
    b2linalg::Vector<T, b2linalg::Vdense> a{index.size()};
    for (size_t i{}; i < index.size(); ++i) { a[i] = dof(index[i], 2); }
 
    f_eff_dyn = m * a;
 
    return std::make_pair(k_eff, f_eff_dyn);
}
 
template <typename T>
typename NewmarkIntegrator<T>::type_t NewmarkIntegrator<T>::type(
      "NewmarkIntegrator", type_name<T>(), StringList(), module, &TimeIntegrator<T>::type);
 
}  // namespace b2000::solver
 
#endif  // B2NEWMARK_H_