b2static_linear_solver.H

//------------------------------------------------------------------------
// b2static_linear_solver.H --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Thomas Blome <thomas.blome@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 B2STATIC_LINEAR_SOLVER_H_
#define B2STATIC_LINEAR_SOLVER_H_
 
#include "b2constraint_util.H"
#include "b2ppconfig.h"
#include "b2solver_system.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 "utils/b2logging.H"
#include "utils/b2object.H"
#include "utils/b2timing.H"
#include "utils/b2util.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class StaticLinearSolver : 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<StaticLinearSolver, Solver::type_t>;
 
    void solve() override {
        CaseList& case_list = model_->get_case_list();
        if (case_list.get_number_of_case() == 1) {
            case_iterator_ = case_list.get_case_iterator();
            solve_on_case(*case_iterator_->next());
        } else if (1) {
            // Solving for each case separately is much
            // faster than solve_cases_optimization().
            case_iterator_ = case_list.get_case_iterator();
            for (;;) {
                auto c = case_iterator_->next();
                if (c == nullptr) { break; }
                solve_on_case(*c);
            }
        }
    }
 
    bool solve_iteration() override {
        solve();
        return false;
    }
 
    static type_t type;
 
protected:
    SolverUtils<T, MATRIX_FORMAT> solver_utile;
    b2linalg::Matrix<T, b2linalg::Mrectangle> f_;  // Contains dof solution. Size: ndofs, nsubcases
    b2linalg::Matrix<T, b2linalg::Mrectangle> f_reac_;
    SolverSystem<T, MATRIX_FORMAT> static_lin_sys;
    void solve_on_case(Case& case_);
};
 
template <typename T, typename MATRIX_FORMAT>
typename StaticLinearSolver<T, MATRIX_FORMAT>::type_t StaticLinearSolver<T, MATRIX_FORMAT>::type(
      "StaticLinearSolver", type_name<T, MATRIX_FORMAT>(), StringList("LINEAR"), module,
      &Solver::type);
 
}  // namespace b2000::solver
 
/* Solve the static linear problems using the Lagrange Multiplier
   Adjunction method.
 
minimize the function E(x) on the variables x (the potential energy
function in a displacement based formulation) defined by:
 
  d_E / d_x (x) = K_i * x + f_i - K_e * x - f_e
 
with the linear constraint on x:
 
  x = R * x_r + r
  L * x + l = 0
 
Where:
  - x (vector) is the degree of freedom (the displacement in a
    displacement based formulation).
  - K_i (matrix) and f_i (vector) define the linear value of the
    elements (the internal force in a displacement based formulation).
  - K_e (matrix) and f_e (vector) define the linear value of the
    natural boundary condition (the external force in a displacement
    based formulation).
  - x_r (vector) is the degree of freedom in the reduction model.
  - R (matrix) and r (vector) define the linear model reduction.
  - L (matrix) and l (vector) define the linear essential boundary conditions.
 
    Using the model reduction equation, the function to minimize become:
 
    trans(R * x_r) * ((K_i - K_e) * R * x_r - f_e + f_i + (K_i - K_e) * r)
 
    with the linear constraint:
 
    L * R * x_r =  -L * r - l
 
    Using the Lagrange Multiplier Adjunction method, the minimization
    problem is equivalent to solve the system of equations:
 
    | trans(R) * (K_i - K_e) * R         trans(L * R) |   | x_r |
    |                                                 | * |     | =
    | L * R                                         0 |   | m   |
 
 
    | trans(R) * (f_e - f_i - (K_i - K_e) * r) |
    |                                          |
    | -L * r - l                               |
 
    The solution x is obtained from the solution x_r using the model
    reduction equation: x = R * x_r + r
 
    The gradients can then be computed.
 
    The difference of the value of the natural boundary condition with
    the value of the element (the reaction force in a displacement
    based formulation) noted f_r is obtained by:
 
    f_r = (K_i * x + f_i) - (K_e * x + f_e)
 
*/
 
template <typename T, typename MATRIX_FORMAT>
void b2000::solver::StaticLinearSolver<T, MATRIX_FORMAT>::solve_on_case(Case& case_) {
    logging::Logger& logger = logging::get_logger("solver.static_linear");
 
    auto time_analysis = obtain_timer("static_linear_analysis");
    auto time_sysresolve = obtain_timer("static_linear_analysis.sysresolve");
 
    (*time_analysis).second.start();
 
    // logger.LOG(logging::info, "Start static linear solver");
    logging::Logger& log_el = logging::get_logger("elementary_matrix");
 
    if (model_ == nullptr) { Exception() << THROW; }
 
    const std::string constraint_string = case_.get_string("CONSTRAINT_METHOD", "REDUCTION");
 
    model_->set_case(case_);
 
    logger << logging::info << "Start static linear solver for case " << case_.get_id()
           << logging::LOGFLUSH;
 
    Domain& domain = model_->get_domain();
    size_t number_of_dof = domain.get_number_of_dof();
    SolutionStaticLinear<T>& solution =
          dynamic_cast<SolutionStaticLinear<T>&>(model_->get_solution());
 
    typename SolutionStaticLinear<T>::gradient_container gc = solution.get_gradient_container();
 
    const bool compute_rcfo =
          solution.compute_constraint_value()
          && (constraint_string == "REDUCTION" || case_.get_bool("GRADIENTS_ONLY", false));
 
    // For the time being no subcases.
    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();
 
    RTable solution_attr;
    f_.resize(number_of_dof, number_of_subcases);
    f_.set_zero();
    if (!case_.get_bool("GRADIENTS_ONLY", false)) {
        static_lin_sys.assemble(model_, case_, solver_utile, solution_attr, f_, f_reac_);
 
        logger.LOG(logging::info, "Resolve the linear problem");
        (*time_sysresolve).second.start();
 
        //
        // Check matrix for structural singularities. Fix them if
        // autospc is active.
        //
        handle_singularities<T, MATRIX_FORMAT>(
              model_, case_, static_lin_sys.K_augm, static_lin_sys.trans_R);
 
        //
        // System resolution.
        //
        const int nbnsv = case_.get_int("COMPUTE_NULL_SPACE_VECTOR", 0);
        b2linalg::Matrix<T, b2linalg::Mrectangle> nks_r;
        static_lin_sys.K_augm.remove_zero();  // Also transforms K_augm to CSC format (TB)
        b2linalg::SingularMatrixError ee;
        bool is_singular = false;
        try {
            static_lin_sys.f_r = inverse(static_lin_sys.K_augm, nks_r, nbnsv) * static_lin_sys.f_r;
        } catch (b2linalg::SingularMatrixError& e) {
            ee = e;
            is_singular = true;
            if (!e.singular_dofs.empty()) {
                if (nks_r.size2() == 0) {
                    nks_r.resize(static_lin_sys.trans_R.size1(), 1);
                    nks_r.set_zero();
                }
                for (auto i : e.singular_dofs) {
                    assert(i < nks_r.size1());
                    for (size_t j = 0; j != nks_r.size2(); ++j) { nks_r(i, j) = T(1); }
                }
            }
        }
        (*time_sysresolve).second.stop();
        solution_attr.set("ELAPSED_TIME_SYSTEM_RESOLUTION", (*time_sysresolve).second.duration());
        if (nks_r.size2() != 0) {
            b2linalg::Matrix<T, b2linalg::Mrectangle> nks;
            nks = transposed(static_lin_sys.trans_R)
                  * nks_r(
                        b2linalg::Interval(0, static_lin_sys.trans_R.size1()),
                        b2linalg::Interval(0, nks_r.size2()));
            solution.set_system_null_space(nks);
        }
        if (is_singular) { ee << THROW; }
 
        // store the dof solution
        f_ = transposed(static_lin_sys.trans_R)
             * static_lin_sys.f_r(
                   b2linalg::Interval(0, static_lin_sys.trans_R.size1()),
                   b2linalg::Interval(0, number_of_subcases));
        for (size_t i = 0; i != number_of_subcases; ++i) {
            f_[i] += static_lin_sys.r;
            solution.set_subcase_id(subcase_ids[i]);
            solution.set_dof(f_[i]);
        }
        solution.set_subcase_id(0);
 
        // computation of the constraint value solution using the solution
        if (solution.compute_constraint_value() && !compute_rcfo) {
            for (size_t i = 0; i != number_of_subcases; ++i) {
                b2linalg::Vector<T, b2linalg::Vdense> ll;
                ll = static_lin_sys.l;
 
                if (constraint_string == "PENALTY" || constraint_string == "AUGMENTED_LAGRANGE") {
                    ll += transposed(static_lin_sys.trans_L) * f_[i];
                    const double penalty_factor = case_.get_double(
                          "CONSTRAINT_PENALTY", constraint_string == "PENALTY" ? 1e10 : 1);
                    ll *= -penalty_factor;
                } else {
                    ll.set_zero();
                }
                if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
                    const double lagrange_mult_scale =
                          case_.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
                    ll += T(-lagrange_mult_scale)
                          * static_lin_sys.f_r[i](b2linalg::Interval(
                                static_lin_sys.trans_R.size1(), static_lin_sys.f_r[i].size()));
                }
                if (case_.get_bool("RCFO_ONLY_SPC", false)) {
                    for (size_t j = 0; j != static_lin_sys.trans_L.size2(); ++j) {
                        if (static_lin_sys.trans_L.get_nb_nonzero(j) != 1) { ll[j] = 0; }
                    }
                }
                f_reac_[i] = static_lin_sys.trans_L * ll;
                solution.set_subcase_id(subcase_ids[i]);
                solution.set_constraint_value(f_reac_[i]);
            }
            solution.set_subcase_id(0);
        }
    } else {
        if (compute_rcfo) {
            // Gradients only: Load the solution field for all subcases.
            for (size_t i = 0; i != subcase_ids.size(); ++i) {
                solution.set_subcase_id(subcase_ids[i]);
                solution.get_dof(f_[i]);
            }
            solution.set_subcase_id(0);
 
            // Get the natural boundary condition (K_nbc_0 * 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>())));
 
            b2linalg::Vector<T, b2linalg::Vdense> f_nbc_0(number_of_dof);
            nbc.get_linear_value(
                  f_nbc_0, b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null);
 
            // Loop over the subcases and add the subcase-specific nbc's.
            for (size_t i = 0; i != subcase_ids.size(); ++i) {
                if (subcase_ids[i] != 0) {
                    nbc.set_subcase_id(subcase_ids[i]);
                    nbc.get_linear_value(
                          f_[i],
                          b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null);
                }
                f_[i] += f_nbc_0;
                solution.set_subcase_id(subcase_ids[i]);
                solution.set_natural_boundary_condition(f_[i]);
            }
            nbc.set_subcase_id(0);
            solution.set_subcase_id(0);
            f_reac_ = -f_;
        }
    }
 
    // Gradients requested: Compute element gradients (stresses, strains,
    // ...) and the constraint value solution ("reaction forces", rcfo)
    if (gc.get() != 0 || compute_rcfo) {
        logger.LOG(logging::info, "Compute gradients and reaction forces");
 
        for (size_t i = 0; i != number_of_subcases; ++i) {
            domain.set_dof(b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(f_[i]));
 
            model_->set_subcase_id(subcase_ids[i]);
 
            Domain::element_iterator ii = domain.get_element_iterator();
            b2linalg::Index index;
            solver_utile.get_elements_values(
                  b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(f_[i]), 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_[i] : 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_[i]); }
        }
        model_->set_subcase_id(0);
 
        (*time_analysis).second.stop();
        solution_attr.set("ELAPSED_TIME_ANALYSIS", (*time_analysis).second.duration());
        solution.terminate_case(true, solution_attr);
    }
 
    logger.LOG(logging::info, "End of static linear solver");
    case_.warn_on_non_used_key();
}
 
#endif  // B2STATIC_LINEAR_SOLVER_H_