b2linearised_prebuckling_solver.H

//------------------------------------------------------------------------
// b2linearised_prebuckling_solver.H --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2004-2012 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (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 B2LINEARIZE_PREBUCKLING_SOLVER_H_
#define B2LINEARIZE_PREBUCKLING_SOLVER_H_
 
#include "b2ppconfig.h"
#include "b2solver_utils.H"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model/b2model.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "solvers/b2constraint_util.H"
#include "utils/b2linear_algebra.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class LinearisedPrebucklingSolver : public Solver {
public:
    using nbc_t = TypedNaturalBoundaryCondition<T, typename MATRIX_FORMAT::sparse_inversible>;
    using ebc_t = TypedEssentialBoundaryCondition<T>;
    using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using type_t = ObjectTypeComplete<LinearisedPrebucklingSolver, Solver::type_t>;
 
    void solve() override {
        while (solve_iteration()) { ; }
    }
 
    bool solve_iteration() override {
        if (model_ == nullptr) { Exception() << THROW; }
        if (case_iterator_.get() == nullptr) {
            case_iterator_ = model_->get_case_list().get_case_iterator();
        }
        case_ = case_iterator_->next();
        if (case_) {
            if (case_->get_number_of_stage() != 1) {
                logging::get_logger("solver.linearised_prebuckling")
                      << logging::warning
                      << "The linearised prebuckling  solver cannot hand a case with multiple "
                         "stages. The extra stages are ignored."
                      << logging::LOGFLUSH;
            }
            solve_on_case(*case_);
            return true;
        }
        return false;
    }
 
    static type_t type;
 
protected:
    void solve_on_case(Case& case_);
    SolverUtils<T, MATRIX_FORMAT> solver_utile;
    // Contains dof solution. Size: ndofs
    b2linalg::Vector<T, b2linalg::Vdense> f_;
    // Number of eigen modes
    int nb_eigenmodes_;
    // Contains eigen value for each mode. Size: nModes
    b2linalg::Vector<T, b2linalg::Vdense> eigenvalues_;
    // Contains eigen vector for each mode. Size: nDofs, nModes
    b2linalg::Matrix<T, b2linalg::Mrectangle> eigenvectors_;
};
 
template <typename T, typename MATRIX_FORMAT>
typename LinearisedPrebucklingSolver<T, MATRIX_FORMAT>::type_t
      LinearisedPrebucklingSolver<T, MATRIX_FORMAT>::type(
            "LinearisedPrebucklingSolver", type_name<T, MATRIX_FORMAT>(),
            StringList("LINEARISED_PREBUCKLING", "L_BIFURCATION"), module, &Solver::type);
 
}  // namespace b2000::solver
 
template <typename T, typename MATRIX_FORMAT>
void b2000::solver::LinearisedPrebucklingSolver<T, MATRIX_FORMAT>::solve_on_case(Case& case_) {
    logging::Logger& logger = logging::get_logger("solver.linearised_prebuckling");
    logging::Logger& log_el = logging::get_logger("elementary_matrix");
 
    if (model_ == nullptr) { Exception() << THROW; }
 
    // Compute the static linear solution
    logger << logging::info << "Start the linearised prebuckling solver for the case "
           << case_.get_id() << "." << logging::LOGFLUSH;
    model_->set_case(case_);
 
    logger << logging::info << "Assemble the linear problem." << logging::LOGFLUSH;
 
    Domain& domain = model_->get_domain();
    size_t number_of_dof = domain.get_number_of_dof();
    SolutionLinearisedPrebuckling<T>& solution =
          dynamic_cast<SolutionLinearisedPrebuckling<T>&>(model_->get_solution());
 
    double sigma;
    char which[3];
    solution.which_modes(nb_eigenmodes_, which, sigma);
    double h = case_.get_double("DLAMBDA", 0.0001);
    const int finite_difference_method = case_.get_int("FINITE_DIFFERENCE_METHOD", 1);
    if (sigma == 0) { sigma = h; }
    const std::string constraint_string = case_.get_string("CONSTRAINT_METHOD", "REDUCTION");
    const bool compute_rcfo =
          solution.compute_constraint_value() && (constraint_string == "REDUCTION");
 
    // Get the model reduction (x = R * x_r + r)
    mrbc_t& mrbc = dynamic_cast<mrbc_t&>(
          model_->get_model_reduction_boundary_condition(mrbc_t::type.get_subtype(
                "TypedModelReductionBoundaryConditionListComponent" + type_name<T>())));
    b2linalg::Vector<T, b2linalg::Vdense> r(number_of_dof);
    b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_R;
    mrbc.get_linear_value(r, trans_R);
 
    // Get the essential boundary condition (L * x + l = 0)
    ebc_t& ebc =
          dynamic_cast<ebc_t&>(model_->get_essential_boundary_condition(ebc_t::type.get_subtype(
                "TypedEssentialBoundaryConditionListComponent" + type_name<T>())));
    b2linalg::Vector<T, b2linalg::Vdense> l(ebc.get_size(true));
    b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_L;
    ebc.get_linear_value(l, trans_L);
 
    if (!l.empty()) {
        // scale the ebc constraint
        b2linalg::Vector<double> trans_L_scale(trans_L.size2());
        trans_L.scale_col_to_norminf(trans_L_scale);
        for (size_t i = 0; i != l.size(); ++i) { l[i] *= trans_L_scale[i]; }
        trans_L.remove_zero();
 
        // procedure to remove dependency on the Lagrange multiplier
        if (case_.get_bool("REMOVE_DEPENDENT_CONSTRAINTS", true)
            && (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE")) {
            b2linalg::Index P;
            // trans_L.get_dep_columns(P, 3e-13, 1.5);
            get_dependent_columns(trans_L, P);
            if (!P.empty()) {
                l.remove(P);
                trans_L.remove_col(P);
                logger << logging::info << "Remove " << P.size() << " dependent constraints."
                       << logging::LOGFLUSH;
            }
        }
    }
 
    // Get the natural boundary condition (K_e * x + f)
    nbc_t& nbc =
          dynamic_cast<nbc_t&>(model_->get_natural_boundary_condition(nbc_t::type.get_subtype(
                "TypedNaturalBoundaryConditionListComponent"
                + type_name<T, typename MATRIX_FORMAT::sparse_inversible>())));
 
    f_.resize(number_of_dof);
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K_e;
    nbc.get_linear_value(f_, K_e);
    solution.set_natural_boundary_condition(f_);
 
    b2linalg::Vector<T, b2linalg::Vdense> f_reac;
    if (compute_rcfo) { f_reac = -f_; }
 
    // Assembly the matrix K by computing
    // (trans(R) * K * R) = K_augm = sum (trans(R) * K_element * R)
    // K * r = K_prod_r = sum (K_element * r)
    logger << logging::info << "Element matrix assembly" << logging::LOGFLUSH;
    size_t K_augm_size = trans_R.size1() + (constraint_string == "PENALTY" ? 0 : trans_L.size2());
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K_augm(
          K_augm_size, domain.get_dof_connectivity_type(), case_);
 
    f_ *= -1;
    solver_utile.get_elements_values(
          b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(r), 1, 0, true, true, trans_R,
          b2linalg::Vector<T, b2linalg::Vdense_constref>::null,
          b2linalg::Vector<double, b2linalg::Vdense_constref>::null, *model_, f_, K_augm,
          b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, log_el);
 
    if (K_e.size1() != 0) {
        K_augm += -trans_R * K_e * transposed(trans_R);
        f_ += -K_e * r;
    }
 
    b2linalg::Vector<T, b2linalg::Vdense> f_r(K_augm_size);
    f_r(b2linalg::Interval(0, trans_R.size1())) = -trans_R * f_;
 
    if (trans_L.size2() > 0) {
        b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
        trans_LR = trans_R * trans_L;
        if (constraint_string == "PENALTY" || constraint_string == "AUGMENTED_LAGRANGE") {
            const double penalty_factor =
                  case_.get_double("CONSTRAINT_PENALTY", constraint_string == "PENALTY" ? 1e10 : 1);
            f_r(b2linalg::Interval(0, trans_R.size1())) += (trans_LR * l) * -penalty_factor;
            K_augm += (trans_LR * transposed(trans_LR)) * penalty_factor;
        }
        if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
            const double lagrange_mult_scale = case_.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
            f_r(b2linalg::Interval(trans_R.size1(), K_augm_size)) = transposed(-trans_L) * r;
            f_r(b2linalg::Interval(trans_R.size1(), K_augm_size)) -= l;
            f_r(b2linalg::Interval(trans_R.size1(), K_augm_size)) *= lagrange_mult_scale;
            trans_LR *= lagrange_mult_scale;
            if (!MATRIX_FORMAT::symmetric) {
                K_augm(
                      b2linalg::Interval(0, trans_R.size1()),
                      b2linalg::Interval(trans_R.size1(), K_augm_size)) += trans_LR;
            }
            b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
            K_augm(
                  b2linalg::Interval(trans_R.size1(), K_augm_size),
                  b2linalg::Interval(0, trans_R.size1())) += LR;
        }
    }
 
    if (logger.is_enabled_for(logging::data)) {
        logger << logging::data << "d_value_d_dof "
               << logging::DBName("D_VALUE_D_DOF.GLOB", 0, 0, case_.get_id()) << K_augm
               << logging::LOGFLUSH;
        logger << logging::data << "value " << logging::DBName("VALUE.GLOB", 0, 0, case_.get_id())
               << f_r << logging::LOGFLUSH;
    }
 
    // system resolution
    logger.LOG(logging::info, "Resolve the linear problem");
    const int nbnsv = case_.get_int("COMPUTE_NULL_SPACE_VECTOR", 0);
    b2linalg::Matrix<T, b2linalg::Mrectangle> nks_r;
    try {
        f_r = inverse(K_augm, nks_r, nbnsv) * f_r;
    } catch (b2linalg::SingularMatrixError& e) {
        if (nks_r.size2() != 0) {
            b2linalg::Matrix<T, b2linalg::Mrectangle> nks;
            nks = transposed(trans_R)
                  * nks_r(
                        b2linalg::Interval(0, trans_R.size1()),
                        b2linalg::Interval(0, nks_r.size2()));
            solution.set_system_null_space(nks);
        }
        throw;
    }
    if (nks_r.size2() != 0) {
        b2linalg::Matrix<T, b2linalg::Mrectangle> nks;
        nks = transposed(trans_R)
              * nks_r(b2linalg::Interval(0, trans_R.size1()), b2linalg::Interval(0, nks_r.size2()));
        solution.set_system_null_space(nks);
    }
 
    // store the dof solution
    f_ = transposed(trans_R) * f_r(b2linalg::Interval(0, trans_R.size1()));
    f_ += r;
    solution.set_dof(f_);
 
    // computation of the constraint value solution using the solution
    if (solution.compute_constraint_value() && !compute_rcfo) {
        if (constraint_string == "PENALTY" || constraint_string == "AUGMENTED_LAGRANGE") {
            l += transposed(trans_L) * f_;
            const double penalty_factor =
                  case_.get_double("CONSTRAINT_PENALTY", constraint_string == "PENALTY" ? 1e10 : 1);
            l *= -penalty_factor;
        } else {
            l.set_zero();
        }
        if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
            const double lagrange_mult_scale = case_.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
            l += T(-lagrange_mult_scale) * f_r(b2linalg::Interval(trans_R.size1(), f_r.size()));
        }
        if (case_.get_bool("RCFO_ONLY_SPC", false)) {
            for (size_t j = 0; j != trans_L.size2(); ++j) {
                if (trans_L.get_nb_nonzero(j) != 1) { l[j] = 0; }
            }
        }
        f_reac = trans_L * l;
        solution.set_constraint_value(f_reac);
    }
 
    typename SolutionStaticLinear<T>::gradient_container gc = solution.get_gradient_container();
 
    // computation of the stress and the constraint value solution
    if (gc.get() != 0 || compute_rcfo) {
        logger.LOG(logging::info, "Compute gradients and reaction forces");
 
        domain.set_dof(b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(f_));
 
        Domain::element_iterator ii = domain.get_element_iterator();
        b2linalg::Index index;
        solver_utile.get_elements_values(
              b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(f_), 1, 0, true, true,
              b2linalg::Matrix<T, b2linalg::Mcompressed_col>::null,
              b2linalg::Vector<T, b2linalg::Vdense_constref>::null,
              b2linalg::Vector<double, b2linalg::Vdense_constref>::null, *model_,
              compute_rcfo ? f_reac : b2linalg::Vector<T, b2linalg::Vdense>::null,
              b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
              b2linalg::Vector<T, b2linalg::Vdense_ref>::null, gc.get(), 0, log_el);
        if (compute_rcfo) { solution.set_constraint_value(f_reac); }
        RTable attributes;
        solution.terminate_case(true, attributes);
    }
 
    logger << logging::info << "Assemble the stability matrix." << logging::LOGFLUSH;
 
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> KG_augm(
          K_augm_size, domain.get_dof_connectivity_type(), case_);
    {
        static const double finite_difference_coefficients[4][4] = {
              {1.0 / 2},
              {2.0 / 3, -1.0 / 12},
              {3.0 / 4, -3.0 / 20, 1.0 / 60},
              {4.0 / 5, -1.0 / 5, 4.0 / 105, -1.0 / 280}};
        if (finite_difference_method < 1 || finite_difference_method > 4) {
            Exception() << "Invalid finite difference method" << THROW;
        }
 
        b2linalg::Vector<double, b2linalg::Vdense> d_value_d_dof_coeff(1);
        for (int i = 0; i != finite_difference_method; ++i) {
            f_ *= i == 0 ? -h : -double(i + 1) / i;
            for (int j = 0; j != 2; j += 1) {
                if (j == 1) { f_ *= -1; }
                d_value_d_dof_coeff[0] =
                      -finite_difference_coefficients[finite_difference_method - 1][i]
                      * (j == 0 ? -1 : 1) / h;
                solver_utile.get_elements_values(
                      b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(f_),
                      h * (i + 1) * (j == 0 ? -1 : 1), 0, true, false, trans_R,
                      b2linalg::Vector<T, b2linalg::Vdense_constref>::null, d_value_d_dof_coeff,
                      *model_, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, KG_augm,
                      b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, log_el);
            }
        }
    }
 
    if (logger.is_enabled_for(logging::data)) {
        logger << logging::data << "d2_value_d_dof_d_lambda "
               << logging::DBName("D2_VALUE_D_DOF_D_LAMBDA.GLOB", 0, 0, case_.get_id()) << KG_augm
               << logging::LOGFLUSH;
    }
 
    if (trans_L.size2() > 0
        && (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE")) {
        // We move the Lagrange multiplier on the KG_augm matrix
        // to have a positive semi definit K_augm matrix
        // for the Buckling spectral shift transformation
        b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
        trans_LR = trans_R * trans_L;
        const double lagrange_mult_scale = case_.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
        trans_LR *= -lagrange_mult_scale;
        if (!MATRIX_FORMAT::symmetric) {
            K_augm(
                  b2linalg::Interval(0, trans_R.size1()),
                  b2linalg::Interval(trans_R.size1(), K_augm_size)) += trans_LR;
        }
        b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
        K_augm(
              b2linalg::Interval(trans_R.size1(), K_augm_size),
              b2linalg::Interval(0, trans_R.size1())) += LR;
        trans_LR *= sigma;
        if (!MATRIX_FORMAT::symmetric) {
            K_augm(
                  b2linalg::Interval(0, trans_R.size1()),
                  b2linalg::Interval(trans_R.size1(), K_augm_size)) += trans_LR;
        }
        LR *= sigma;
        K_augm(
              b2linalg::Interval(trans_R.size1(), K_augm_size),
              b2linalg::Interval(0, trans_R.size1())) += LR;
    }
 
    logger << logging::info << "Eigenvalue problem resolution" << logging::LOGFLUSH;
 
    eigenvalues_.resize(nb_eigenmodes_);
    b2linalg::Matrix<T, b2linalg::Mrectangle> eigenvector_red(K_augm_size, nb_eigenmodes_);
 
    eigenvector_red = eigenvector(K_augm, 'P', KG_augm, 'H', eigenvalues_, which, sigma);
 
    eigenvectors_.resize(number_of_dof, nb_eigenmodes_);
    mrbc.get_linear_value(b2linalg::Vector<T, b2linalg::Vdense>::null, trans_R);
    eigenvectors_ =
          transposed(trans_R)
          * eigenvector_red(
                b2linalg::Interval(0, trans_R.size1()), b2linalg::Interval(0, nb_eigenmodes_));
 
    solution.set_modes(eigenvalues_, eigenvectors_);
    RTable attributes;
    solution.terminate_case(true, attributes);
    case_.warn_on_non_used_key();
    logger.LOG(logging::info, "End of linearised prebuckling solver");
}
 
#endif  // B2LINEARIZE_PREBUCKLING_SOLVER_H_