b2solver_system.H

//------------------------------------------------------------------------
// b2solver_system.H --
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2004-2012,2016,2017 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 B2SOLVER_SYSTEM_H_
#define B2SOLVER_SYSTEM_H_
 
#include "b2constraint_util.H"
#include "b2ppconfig.h"
#include "b2solver_utils.H"
#include "model/b2boundary_condition.H"
#include "model/b2case.H"
#include "model/b2domain.H"
#include "model/b2model.H"
#include "model/b2solution.H"
#include "utils/b2logging.H"
#include "utils/b2timing.H"
#include "utils/b2util.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class SolverSystem {
public:
    void assemble(
          Model* const model, Case const& case_, SolverUtils<T, MATRIX_FORMAT>& solver_utile,
          RTable& solution_attr, b2linalg::Matrix<T, b2linalg::Mrectangle>& f_,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& f_reac);
 
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K_augm;
    b2linalg::Matrix<T, b2linalg::Mrectangle> f_r;
    b2linalg::Vector<T, b2linalg::Vdense> l;
    b2linalg::Vector<T, b2linalg::Vdense> r;
    b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_L;
    b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_R;
    b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense> f_nbc_0;
};
 
template <typename T>
void bc_elimination(
      Model* const model, Case const& case_, b2linalg::Vector<T, b2linalg::Vdense>& r,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_R,
      b2linalg::Vector<T, b2linalg::Vdense>& l,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_L) {
    using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using ebc_t = TypedEssentialBoundaryCondition<T>;
 
    logging::Logger& logger = logging::get_logger("solver.static_linear");
 
    size_t number_of_dof = model->get_domain().get_number_of_dof();
 
    // 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>())));
    r = b2linalg::Vector<T, b2linalg::Vdense>(number_of_dof);
    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>())));
    l = b2linalg::Vector<T, b2linalg::Vdense>(ebc.get_size(true));
    ebc.get_linear_value(l, trans_L);
 
    if (!l.empty()) {
        const std::string constraint_string = case_.get_string("CONSTRAINT_METHOD", "REDUCTION");
        // 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 << "Removed " << P.size() << " dependent constraints."
                       << logging::LOGFLUSH;
            }
        }
    }
 
}  // bc_elimination
 
template <typename T, typename MATRIX_FORMAT>
void apply_nbcs(
      Model* const model, b2linalg::Matrix<T, b2linalg::Mrectangle>& f_,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& f_reac,
      b2linalg::Vector<T, b2linalg::Vdense>& f_nbc_0,
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_nbc_0) {
    using nbc_t = TypedNaturalBoundaryCondition<T, typename MATRIX_FORMAT::sparse_inversible>;
 
    logging::Logger& logger = logging::get_logger("solver.static_linear");
 
    std::set<int> subcase_id_set;
    model->get_subcase_ids(subcase_id_set);
    std::vector<int> subcase_ids(subcase_id_set.begin(), subcase_id_set.end());
    size_t number_of_dof = model->get_domain().get_number_of_dof();
    SolutionStaticLinear<T>& solution =
          dynamic_cast<SolutionStaticLinear<T>&>(model->get_solution());
 
    logger << logging::info << "Assemble the natural boundary conditions" << logging::LOGFLUSH;
    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>())));
 
    // Calculate the natural bc's that are valid for all subcases.
    f_nbc_0 = b2linalg::Vector<T, b2linalg::Vdense>(number_of_dof);
    nbc.get_linear_value(f_nbc_0, K_nbc_0);
 
    // Loop over the subcases and add the subcase-specific nbc's to
    // the vector of forces.
    for (size_t i = 0; i != subcase_ids.size(); ++i) {
        model->set_subcase_id(subcase_ids[i]);
        if (subcase_ids[i] != 0) {
            // Missing: add check whether nbc gives a second variation,
            // in this case, an exception must be thrown.
            nbc.get_linear_value(
                  f_[i], b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null);
        }
        f_[i] += f_nbc_0;
        solution.set_natural_boundary_condition(f_[i]);
    }
    model->set_subcase_id(0);
    f_reac = -f_;
    f_ *= -1;
}  // apply_nbcs
 
template <typename T, typename MATRIX_FORMAT>
void assemble_from_elems(
      Model* const model, SolverUtils<T, MATRIX_FORMAT>& solver_utile,
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> const& K_nbc_0,
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& f_, b2linalg::Matrix<T, b2linalg::Mrectangle>& f_r,
      b2linalg::Vector<T, b2linalg::Vdense>& r,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_R) {
    logging::Logger& log_el = logging::get_logger("elementary_matrix");
 
    Domain& domain = model->get_domain();
    std::set<int> subcase_id_set;
    model->get_subcase_ids(subcase_id_set);
    std::vector<int> subcase_ids(subcase_id_set.begin(), subcase_id_set.end());
    const size_t number_of_subcases = subcase_ids.size();
    size_t number_of_dof = model->get_domain().get_number_of_dof();
 
    b2linalg::Vector<T, b2linalg::Vdense> f_add_to_all(number_of_dof);
    f_add_to_all.set_zero();
    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_add_to_all, K_augm,
          b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, log_el);
 
    if (K_nbc_0.size1() != 0) {
        K_augm += -trans_R * K_nbc_0 * transposed(trans_R);
        f_add_to_all += -K_nbc_0 * r;
    }
    for (size_t i = 0; i != subcase_ids.size(); ++i) {
        model->set_subcase_id(subcase_ids[i]);
        if (subcase_ids[i] != 0 && domain.have_temperature()) {
            // Calculate the element 1st variation for this subcase.
            b2linalg::Vector<T, b2linalg::Vdense> f_add_to_subcase(number_of_dof);
            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_add_to_subcase,
                  b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
                  b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, log_el);
            if (K_nbc_0.size1() != 0) { f_add_to_subcase += -K_nbc_0 * r; }
            f_[i] += f_add_to_subcase;
        } else {
            f_[i] += f_add_to_all;
        }
    }
    model->set_subcase_id(0);
    f_r(b2linalg::Interval(0, trans_R.size1()), b2linalg::Interval(0, number_of_subcases)) =
          -trans_R * f_;
}  // assemble_from_elems
 
template <typename T, typename MATRIX_FORMAT>
void apply_constraints(
      Model* const model, Case const& case_,
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& f_r, b2linalg::Vector<T, b2linalg::Vdense>& l,
      b2linalg::Vector<T, b2linalg::Vdense>& r,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_L,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_R) {
    const std::string constraint_string = case_.get_string("CONSTRAINT_METHOD", "REDUCTION");
    size_t K_augm_size = trans_R.size1() + (constraint_string == "PENALTY" ? 0 : trans_L.size2());
    std::set<int> subcase_id_set;
    model->get_subcase_ids(subcase_id_set);
    std::vector<int> subcase_ids(subcase_id_set.begin(), subcase_id_set.end());
    const size_t number_of_subcases = subcase_ids.size();
 
    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);
            for (size_t i = 0; i != number_of_subcases; ++i) {
                f_r[i](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);
            for (size_t i = 0; i != number_of_subcases; ++i) {
                f_r[i](b2linalg::Interval(trans_R.size1(), K_augm_size)) = transposed(-trans_L) * r;
                f_r[i](b2linalg::Interval(trans_R.size1(), K_augm_size)) -= l;
            }
            f_r(b2linalg::Interval(trans_R.size1(), K_augm_size),
                b2linalg::Interval(0, number_of_subcases)) *= 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;
        }
    }
}  // apply_constraints
 
template <typename T, typename MATRIX_FORMAT>
void handle_singularities(
      Model* const model, Case const& case_,
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_R) {
    SolutionStaticLinear<T>& solution =
          dynamic_cast<SolutionStaticLinear<T>&>(model->get_solution());
    const bool autospc = case_.get_bool("AUTOSPC", false);
 
    std::set<size_t> zero_columns_r;
    for (size_t i = 0; i != K_augm.size1(); ++i) {
        if (K_augm(i, i) == T(0)) {
            bool all_zero = true;
            for (size_t j = 0; all_zero && j != K_augm.size2(); ++j) {
                if (K_augm(i, j) != T(0)) { all_zero = false; }
            }
            if (all_zero) {
                zero_columns_r.insert(i);
                if (autospc) { K_augm(i, i) = 1.; }
            }
        }
    }
 
    if (!zero_columns_r.empty()) {
        b2linalg::Matrix<T, b2linalg::Mrectangle> nks_r;
        nks_r.resize(trans_R.size1(), 1);
        nks_r.set_zero();
        for (auto i : zero_columns_r) { nks_r(i, 0) = T(1); }
        b2linalg::Matrix<T, b2linalg::Mrectangle> nks;
        nks = transposed(trans_R)
              * nks_r(b2linalg::Interval(0, trans_R.size1()), b2linalg::Interval(0, nks_r.size2()));
        for (size_t i = 0; i != nks.size1(); ++i) { nks(i, 0) = (nks(i, 0) != T(0) ? T(1) : T(0)); }
        solution.set_system_null_space(nks, false, "NS_STRUCTURE");
    }
}  // handle_singularities
 
// Assemble the system for a static linear problem
template <typename T, typename MATRIX_FORMAT>
void assemble_system(
      Model* const model, Case const& case_, SolverUtils<T, MATRIX_FORMAT>& solver_utile,
      RTable& solution_attr, b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& f_, b2linalg::Matrix<T, b2linalg::Mrectangle>& f_r,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& f_reac, b2linalg::Vector<T, b2linalg::Vdense>& l,
      b2linalg::Vector<T, b2linalg::Vdense>& r, b2linalg::Vector<T, b2linalg::Vdense>& f_nbc_0,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_L,
      b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_R) {
    logging::Logger& logger = logging::get_logger("solver.static_linear");
    auto time_nbc = obtain_timer("static_linear_analysis.nbc");
    auto time_elements = obtain_timer("static_linear_analysis.elements");
 
    const std::string constraint_string = case_.get_string("CONSTRAINT_METHOD", "REDUCTION");
    Domain& domain = model->get_domain();
    std::set<int> subcase_id_set;
    model->get_subcase_ids(subcase_id_set);
    std::vector<int> subcase_ids(subcase_id_set.begin(), subcase_id_set.end());
    const size_t number_of_subcases = subcase_ids.size();
 
    bc_elimination<T>(model, case_, r, trans_R, l, trans_L);
 
    // 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)
 
    // Apply natural boundary conditions
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K_nbc_0;
    (*time_nbc).second.start();
    apply_nbcs<T, MATRIX_FORMAT>(model, f_, f_reac, f_nbc_0, K_nbc_0);
    (*time_nbc).second.stop();
    solution_attr.set("ELAPSED_TIME_NBC", (*time_nbc).second.duration());
 
    // Add the element contribution and the dof-dependent nbc.
    logger << logging::info << "Element matrix assembly" << logging::LOGFLUSH;
    (*time_elements).second.start();
 
    size_t K_augm_size = trans_R.size1() + (constraint_string == "PENALTY" ? 0 : trans_L.size2());
    K_augm = b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>(
          K_augm_size, domain.get_dof_connectivity_type(), case_);
    f_r = b2linalg::Matrix<T, b2linalg::Mrectangle>(K_augm_size, number_of_subcases);
    assemble_from_elems<T, MATRIX_FORMAT>(
          model, solver_utile, K_nbc_0, K_augm, f_, f_r, r, trans_R);
 
    // Apply the constraints to the system
    apply_constraints<T, MATRIX_FORMAT>(model, case_, K_augm, f_r, l, r, trans_L, trans_R);
 
    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;
    }
    (*time_elements).second.stop();
    solution_attr.set("ELAPSED_TIME_ELEMENTS", (*time_elements).second.duration());
 
}  // assemble_system
 
template <typename T, typename MATRIX_FORMAT>
void SolverSystem<T, MATRIX_FORMAT>::assemble(
      Model* const model, Case const& case_, SolverUtils<T, MATRIX_FORMAT>& solver_utile,
      RTable& solution_attr, b2linalg::Matrix<T, b2linalg::Mrectangle>& f_,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& f_reac) {
    assemble_system<T, MATRIX_FORMAT>(
          model, case_, solver_utile, solution_attr, K_augm, f_, f_r, f_reac, l, r, f_nbc_0,
          trans_L, trans_R);
}
 
}  // namespace b2000::solver
 
#endif  // B2SOLVER_SYSTEM_H_