b2boundary_condition_nlr_coupling.H

//------------------------------------------------------------------------
// b2boundary_condition_database.C --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 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 B2BOUNDARY_CONDITION_NLR_COUPLING_H_
#define B2BOUNDARY_CONDITION_NLR_COUPLING_H_
 
#include <set>
 
#include "b2model_database.H"
#include "b2solution_database.H"
#include "model/b2case.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model_imp/b2boundary_condition_component.H"
#include "utils/b2linear_algebra.H"
 
// #define BCNC_CHECK_DERIVATIVE 1
// #define CHECK_GLOBAL_TO_LOCAL_MAPPING 1e-5
 
namespace b2000::b2dbv3 {
 
class StaticResidueFunctionNLRCoupling {
public:
    typedef double T;
    typedef b2linalg::Msymmetric MATRIX_FORMAT;
    typedef TypedNaturalBoundaryConditionListComponent<T, MATRIX_FORMAT::sparse_inversible>
          nbc_list_t;
    typedef TypedEssentialBoundaryConditionListComponent<T> ebc_list_t;
    typedef TypedModelReductionBoundaryConditionListComponent<T> mrbc_list_t;
 
    StaticResidueFunctionNLRCoupling()
        : number_of_dof(0),
          number_of_non_lagrang_dof(0),
          model(nullptr),
          domain(nullptr),
          mrbc(nullptr),
          ebc(nullptr),
          nbc(nullptr),
          constraint_has_penalty(false),
          penalty_factor(0),
          lagrange_mult_scale(0) {}
 
    virtual ~StaticResidueFunctionNLRCoupling() {}
 
    void init(Model& model_, const AttributeList& attribute) {
        model = &model_;
        domain = &model->get_domain();
        ebc = &dynamic_cast<ebc_list_t::typed_base_t&>(
              model->get_essential_boundary_condition(&ebc_list_t::type));
        mrbc = &dynamic_cast<mrbc_list_t::typed_base_t&>(
              model->get_model_reduction_boundary_condition(&mrbc_list_t::type));
        nbc = &dynamic_cast<nbc_list_t::typed_base_t&>(
              model->get_natural_boundary_condition(&nbc_list_t::type));
        nbc_sol.resize(domain->get_number_of_dof());
        fi.resize(domain->get_number_of_dof());
 
        constraint_has_penalty =
              attribute.get_string("CONSTRAINT_METHOD", "REDUCTION") == "PENALTY";
 
        constraint_has_penalty |=
              attribute.get_string("NONLINEAR_CONSTRAINT_METHOD", "OTHER") == "PENALTY";
 
        number_of_dof = number_of_non_lagrang_dof =
              mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
        if (constraint_has_penalty) {
            penalty_factor = attribute.get_double("CONSTRAINT_PENALTY", 1e10);
        } else {
            number_of_dof += ebc->get_size(false);
            penalty_factor = attribute.get_double("CONSTRAINT_PENALTY", 1);
            if (attribute.get_string("NONLINEAR_CONSTRAINT_METHOD", "OTHER") == "LAGRANGE") {
                penalty_factor = 0;
            }
            lagrange_mult_scale = attribute.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
        }
    }
 
    size_t get_number_of_dof() { return number_of_dof; }
 
    size_t get_number_of_non_lagrange_dof() { return number_of_non_lagrang_dof; }
 
    void get_solution(
          const b2linalg::Vector<T, b2linalg::Vdense_constref> dof, const double time,
          EquilibriumSolution equilibrium_solution,
          b2linalg::Vector<T, b2linalg::Vdense_ref> dof_sol) {
        mrbc_get_nonlinear_value(
              dof(b2linalg::Interval(
                    0, mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1))),
              time, equilibrium_solution, dof_sol,
              b2linalg::Matrix<T, b2linalg::Mcompressed_col>::null,
              b2linalg::Vector<T, b2linalg::Vdense_ref>::null, nullptr);
    }
 
    const b2linalg::Vector<T, b2linalg::Vdense_constref> get_residuum_sol() { return fi; }
 
    const b2linalg::Vector<T, b2linalg::Vdense_constref> get_nbc_sol() { return nbc_sol; }
 
    T get_residue_normalization() { return std::max(fi * fi, nbc_sol * nbc_sol); }
 
    double get_energy_normalisation1() {
        double res = 0;
        for (size_t i = 0; i != r.size(); ++i) {
            res += std::max(norm(r[i] * nbc_sol[i]), norm(r[i] * fi[i]));
        }
        return res;
    }
 
    double get_energy() {
        double res = 0;
        for (size_t i = 0; i != r.size(); ++i) { res += r[i] * (fi[i] + nbc_sol[i]); }
        return res;
    }
 
    void get_residue(
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof_global,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, const double time,
          const double delta_time, const EquilibriumSolution equilibrium_solution,
          b2linalg::Vector<T, b2linalg::Vdense_ref> residue,
          b2linalg::Matrix<T, MATRIX_FORMAT::sparse_inversible>& d_residue_d_dof,
          const b2linalg::Matrix<T>& C, b2linalg::Matrix<T>& RtKC,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& CtKC,
          b2linalg::Vector<T, b2linalg::Vdense_ref> d_residue_d_time, SolverHints* solver_hints) {
        b2linalg::Interval disp_range(
              0, mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
        b2linalg::Interval lagrange_mult_range(
              mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1), number_of_dof);
 
        const size_t nb_dof = domain->get_number_of_dof();
 
        if (!d_residue_d_dof.is_null()) {
            if (d_residue_d_dof.size1() != dof.size()) { d_residue_d_dof.resize(dof.size()); }
            d_residue_d_dof.set_zero_same_struct();
            RtKC.resize(dof.size(), C.size2());
            RtKC.set_zero();
        }
 
        r.resize(nb_dof);
        b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
        b2linalg::Vector<T, b2linalg::Vdense> d_r_d_time(domain->get_number_of_dof());
        mrbc_get_nonlinear_value(
              dof(disp_range), time, equilibrium_solution, r, trans_d_r_d_dof, d_r_d_time,
              solver_hints);
        r += dof_global;
 
        b2linalg::Matrix<T, MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
        b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(domain->get_number_of_dof());
        nbc_sol.set_zero();
        nbc->get_nonlinear_value(
              r, time, equilibrium_solution, residue,
              (d_residue_d_dof.is_null()
                     ? b2linalg::Matrix<T, MATRIX_FORMAT::sparse_inversible>::null
                     : d_fe_d_dof),
              (d_residue_d_time.is_null() ? d_residue_d_time
                                          : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time)),
              solver_hints);
        fi = -nbc_sol;
        d_fe_d_time = -d_fe_d_time;
        if (d_fe_d_dof.size1() != 0) {
            d_residue_d_dof += -trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof);
        }
 
        b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
        b2linalg::Vector<T> l;
 
        if (constraint_has_penalty) {
            l.resize(ebc->get_size(false));
            ebc->get_nonlinear_value(
                  r, time, equilibrium_solution, l, trans_d_l_d_dof,
                  b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
        } else {
            ebc->get_nonlinear_value(
                  r, time, equilibrium_solution, residue(lagrange_mult_range), trans_d_l_d_dof,
                  d_residue_d_time.is_null() ? d_residue_d_time
                                             : d_residue_d_time(lagrange_mult_range),
                  solver_hints);
            if (!d_residue_d_time.is_null()) {
                d_residue_d_time(lagrange_mult_range) +=
                      b2linalg::transposed(trans_d_l_d_dof) * d_r_d_time;
            }
        }
 
        if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
            d_fe_d_time -= d_fe_d_dof * d_r_d_time;
        }
 
        // iterate on the elements
        {
            domain->set_dof(r);
            b2linalg::Index index;
            b2linalg::Vector<T, b2linalg::Vdense> value_e;
            b2linalg::Vector<T, b2linalg::Vdense> d_value_d_time_e;
            std::vector<bool> d_value_d_dof_e_flags(
                  d_residue_d_dof.is_null() && d_residue_d_time.is_null() ? 0 : 1, true);
            std::vector<b2linalg::Matrix<T, MATRIX_FORMAT::dense> > d_value_d_dof_e(
                  d_value_d_dof_e_flags.size());
            std::unique_ptr<b2000::Domain::ElementIterator> i = domain->get_element_iterator();
            if (K.size1() == 0) { K.resize(nb_dof, nb_dof); }
            K.set_zero_same_struct();
            for (Element* element = i->next(); element != nullptr; element = i->next()) {
                if (auto te = dynamic_cast<TypedElement<T>*>(element); te != nullptr) {
                    te->get_value(
                          *model, false, equilibrium_solution, time, delta_time,
                          b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(r),
                          nullptr,  // gradient_container
                          solver_hints, index, value_e, d_value_d_dof_e_flags, d_value_d_dof_e,
                          (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense>::null
                                                      : d_value_d_time_e));
                }
                fi(index) += value_e;
 
                if (!d_residue_d_time.is_null()) {
                    d_value_d_time_e += d_value_d_dof_e[0] * d_r_d_time(index);
                    d_residue_d_time(index) += d_value_d_time_e;
                }
                if (!d_residue_d_dof.is_null()) {
                    b2linalg::Matrix<
                          T, b2linalg::Mstructured_constref<b2linalg::Mpacked_st_constref> >
                          Kelem = b2linalg::StructuredMatrix(nb_dof, index, d_value_d_dof_e[0]);
                    d_residue_d_dof +=
                          trans_d_r_d_dof * Kelem * b2linalg::transposed(trans_d_r_d_dof);
                    K += Kelem;
                }
            }
            {
                b2linalg::Vector<double> KCi;
                b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref> Kref(K);
                for (size_t j = 0; j != C.size2(); ++j) {
                    KCi = Kref * C[j];
                    RtKC[j] += trans_d_r_d_dof * KCi;
                    CtKC[j] += b2linalg::transposed(C) * KCi;
                }
            }
        }
 
        {
            b2linalg::Vector<T> tmp = fi;
            if (constraint_has_penalty) {
                tmp += (trans_d_l_d_dof * l) * penalty_factor;
            } else {
                if (MATRIX_FORMAT::symmetric && penalty_factor != 0) {
                    // We apply penalty method with a small factor to obtain a definite positive
                    // matrix.
                    b2linalg::Vector<T> tmp1;
                    tmp1 = residue(lagrange_mult_range) * penalty_factor;
                    tmp1 += dof(lagrange_mult_range) * lagrange_mult_scale;
                    tmp += trans_d_l_d_dof * tmp1;
                    residue(lagrange_mult_range) *= lagrange_mult_scale;
                } else {
                    tmp += lagrange_mult_scale * trans_d_l_d_dof * dof(lagrange_mult_range);
                }
            }
            residue(disp_range) = trans_d_r_d_dof * tmp;
        }
 
        if (!d_residue_d_dof.is_null()) {
            b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
            trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
            if (penalty_factor != 0 && (constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
                // We apply penalty method but with a small factor to
                // obtain a definite positive matrix when Lagrange multipliers method is used.
                d_residue_d_dof(disp_range, disp_range) +=
                      (trans_LR * transposed(trans_LR)) * penalty_factor;
            }
            if (!constraint_has_penalty) {
                trans_LR *= lagrange_mult_scale;
                b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
                d_residue_d_dof(lagrange_mult_range, disp_range) += LR;
                if (!MATRIX_FORMAT::symmetric) {
                    d_residue_d_dof(disp_range, lagrange_mult_range) += trans_LR;
                }
            }
        }
 
        if (!d_residue_d_time.is_null()) {
            d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
        }
    }
 
    void mrbc_get_nonlinear_inverse_value(
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, double time,
          b2linalg::Vector<T, b2linalg::Vdense_ref> reduced_dof) {
        mrbc->get_nonlinear_inverse_value(
              b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(dof), time,
              b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(reduced_dof));
    }
 
protected:
    virtual void mrbc_get_nonlinear_value(
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, double time,
          EquilibriumSolution equilibrium_solution, b2linalg::Vector<T, b2linalg::Vdense_ref> value,
          b2linalg::Matrix<T, b2linalg::Mcompressed_col>& d_value_d_dof_trans,
          b2linalg::Vector<T, b2linalg::Vdense_ref> d_value_d_time, SolverHints* solver_hints) {
        mrbc->get_nonlinear_value(
              dof, time, equilibrium_solution, value, d_value_d_dof_trans,
              b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(d_value_d_time), solver_hints);
    }
 
    size_t number_of_dof;
    size_t number_of_non_lagrang_dof;
    Model* model;
    b2000::Domain* domain;
    TypedModelReductionBoundaryCondition<T>* mrbc;
    ebc_list_t::typed_base_t* ebc;
    nbc_list_t::typed_base_t* nbc;
    b2linalg::Vector<T, b2linalg::Vdense> fi;
    b2linalg::Vector<T, b2linalg::Vdense> nbc_sol;
    b2linalg::Vector<T, b2linalg::Vdense> r;
    bool constraint_has_penalty;
    double penalty_factor;
    double lagrange_mult_scale;
    b2linalg::Matrix<T, MATRIX_FORMAT::sparse_inversible> K;
};
 
class NBC_NLR_Coupling;
 
class SolutionStaticNonlinearNLRCoupling : public b2000::SolutionStaticNonlinear<double> {
public:
    SolutionStaticNonlinearNLRCoupling()
        : nlr_coupling(nullptr), local_solution(nullptr), subcase_id(0) {}
 
    void set_nlr_coupling(NBC_NLR_Coupling& nlr_coupling);
 
    void set_model(b2000::Model& model) override {}
 
    void set_case(b2000::Case& case_) override {}
 
    void set_subcase_id(const int subcase_id_) override { subcase_id = subcase_id_; }
 
    void set_need_refinement() override {}
 
    b2000::SolutionStaticNonlinear<double>::gradient_container get_gradient_container() override {
        return local_solution->get_gradient_container();
    }
 
    void new_cycle(bool end_of_stage = false) override { local_solution->new_cycle(end_of_stage); }
 
    void terminate_cycle(double time, double dtime, const RTable& attributes) override {
        local_solution->terminate_cycle(time, dtime, attributes);
    }
 
    void terminate_stage(bool success = true) override;
 
    void terminate_case(bool success, const RTable& attributes) override;
 
    int get_subcase_id() const { return subcase_id; }
 
    size_t get_cycle_id() override { return local_solution->get_cycle_id(); }
 
    size_t get_subcycle_id() override { return local_solution->get_subcycle_id(); }
 
    std::string get_domain_state_id() override { return local_solution->get_domain_state_id(); }
 
    void set_dof(
          double time, const b2linalg::Vector<double, b2linalg::Vdense_constref>& dof,
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& nbc,
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& residue,
          bool unconverged_subcycle = false) override;
 
    void set_equilibrium(
          double time, const b2linalg::Vector<double, b2linalg::Vdense_constref>& dof) override {}
 
    void get_dof(double& time, b2linalg::Vector<double, b2linalg::Vdense_ref> dof) override {
        UnimplementedError() << THROW;
    }
 
private:
    NBC_NLR_Coupling* nlr_coupling;
    b2000::b2dbv3::SolutionStaticNonlinear<double>* local_solution;
    int subcase_id;
};
 
class NBC_NLR_Coupling : public TypedNaturalBoundaryConditionComponent<
                               double, b2linalg::Msym_compressed_col_update_inv> {
public:
    NBC_NLR_Coupling() : local_case(nullptr), global_model(nullptr), C_linear(true) {}
 
    void init(
          const std::string& id_, TimeFunction* scalef_, b2000::Model& model_,
          const b2000::Case& case_, const int subcase_id_) override {
        id = id_;
        global_model = &model_;
        subcase_id = subcase_id_;
    }
 
    std::string get_id() const override { return id; }
 
    int get_subcase_id() const override { return subcase_id; }
 
    NaturalBoundaryCondition::Type get_nbc_type() const override {
        return NaturalBoundaryCondition::conservative;
    }
 
    void add_nonlinear_nbc_value(
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& global_dof_all, double time,
          EquilibriumSolution equilibrium_solution,
          b2linalg::Vector<double, b2linalg::Vdense_ref> value_,
          b2linalg::Matrix<double, b2linalg::Msym_compressed_col_update_inv>& d_value_d_dof,
          b2linalg::Vector<double, b2linalg::Vdense_ref> d_value_d_time,
          SolverHints* solver_hints) override {
        // warning: the ebv values are subtracted to the internal force of the model
 
        if (C.size1() == 0) { update_local_model(); }
 
        if (equilibrium_solution.input_level == 0) {
            prev_global_dof = conv_global_dof;
            prev_local_dof = conv_local_dof;
        }
        b2linalg::Vector<double> local_dof = prev_local_dof;
 
        b2linalg::Vector<double> global_dof = global_dof_all(C_index);
        // update the local dof using the delta of the global dof
        {
            if (prev_global_dof.size() != 0) {
                b2linalg::Vector<double> delta_global_dof;
                delta_global_dof = prev_global_dof - global_dof;
                b2linalg::Vector<double> delta_d;
                delta_d = local_RtKC * delta_global_dof;
                delta_d = b2linalg::inverse(local_RtKR) * delta_d;
                local_dof += delta_d;
            }
            if (equilibrium_solution.output_level == 0) { conv_global_dof = global_dof; }
            if (equilibrium_solution.output_level != -1) { prev_global_dof = global_dof; }
        }
 
        b2linalg::Matrix<double, b2linalg::Mrectangle> local_CtKC(C.size2(), C.size2());
 
        // compute the residue of the local model and this derivative
        b2linalg::Vector<double> residue(local_model_residue.get_number_of_dof());
        {
            b2linalg::Vector<double> global_dof_on_local;
            if (C_linear) {
                global_dof_on_local = C * global_dof;
            } else {
                get_global_to_local_dof(global_dof_all, global_dof_on_local, C);
            }
            local_model_residue.get_residue(
                  global_dof_on_local, local_dof, time, 0, equilibrium_solution, residue,
                  (d_value_d_dof.is_null()
                         ? b2linalg::Matrix<double, b2linalg::Msym_compressed_col_update_inv>::null
                         : local_RtKR),
                  C, local_RtKC, local_CtKC, b2linalg::Vector<double, b2linalg::Vdense_ref>::null,
                  solver_hints);
        }
 
        // compute the contribution to the second variation of the global model
        if (!d_value_d_dof.is_null() && local_RtKR.size1() != 0) {
            b2linalg::Matrix<double> M;
            M = inverse(local_RtKR) * local_RtKC;
            local_CtKC -= (b2linalg::transposed(local_RtKC) * M);
 
            // addition of the second derivative component
            // This is expensive (numerical computation of the derivative and don't improve the
            // convergence
            /*if (!C_linear) {
                b2linalg::Vector<double> global_dof_all_tmp = global_dof_all;
                b2linalg::Matrix<double> C_tmp(C.size1(), C.size2());
                b2linalg::Vector<double> global_dof_on_local_tmp;
                const double h = 1e-3;
                for (size_t i = 0; i != C_index.size(); ++i) {
                    const double dtmp = global_dof_all_tmp[i];
                    global_dof_all_tmp[i] += h;
                    get_global_to_local_dof(global_dof_all_tmp, global_dof_on_local_tmp, C_tmp);
                    C_tmp += -C;
                    {
                        std::cout << i << std::endl;
                        b2linalg::Vector<double> tmp;
                        tmp = ((1 /h) * (transposed(C_tmp) *
            local_model_residue.get_residuum_sol())); std::cout << tmp << std::endl;
                    }
                    local_CtKC[i] += (1 /h) * (transposed(C_tmp) *
            local_model_residue.get_residuum_sol()); global_dof_all_tmp[i] = dtmp;
                }
            }*/
        }
 
        b2linalg::Vector<double> delta_d;
        // update the local dof using the residue (for the next Newton iteration)
        if (equilibrium_solution.output_level != -1 || !value_.is_null()) {
            delta_d = b2linalg::inverse(local_RtKR) * residue;
            if (equilibrium_solution.output_level != -1) {
                local_dof -= delta_d;
                prev_local_dof = local_dof;
            }
        }
 
        // compute the contribution to the first variation of the global model
        if (!value_.is_null()) {
            b2linalg::Vector<double> tmp;
            tmp = -transposed(C) * local_model_residue.get_residuum_sol();
            tmp += transposed(local_RtKC) * delta_d;
            value_(C_index) += tmp;
        }
 
        // compute the contribution to the second variation of the global model
        if (!d_value_d_dof.is_null()) {
            if (d_value_d_dof.size1() == 0) {
                d_value_d_dof.resize(global_dof_all.size(), global_dof_all.size());
            }
            b2linalg::Matrix<double, b2linalg::Mupper_packed> local_CtKC_sym;
            local_CtKC_sym = -local_CtKC;
            d_value_d_dof +=
                  b2linalg::StructuredMatrix(d_value_d_dof.size1(), C_index, local_CtKC_sym);
        }
 
        // remove of the global element contribution
        {
            b2linalg::Index index_e;
            b2linalg::Vector<double, b2linalg::Vdense> value_e;
            b2linalg::Vector<double, b2linalg::Vdense> d_value_d_time_e;
            std::vector<bool> d_value_d_dof_e_flags(
                  d_value_d_dof.is_null() && d_value_d_time.is_null() ? 0 : 1, true);
            std::vector<b2linalg::Matrix<double, b2linalg::Mpacked> > d_value_d_dof_e(
                  d_value_d_dof_e_flags.size());
            for (global_elements_t::iterator e = global_elements.begin();
                 e != global_elements.end(); ++e) {
                if (auto te = dynamic_cast<TypedElement<double>*>(*e); te != nullptr) {
                    te->get_value(
                          *global_model, false, equilibrium_solution, time, 0,
                          b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(global_dof_all),
                          nullptr,  // gradient_container
                          solver_hints, index_e,
                          (value_.is_null() ? b2linalg::Vector<double, b2linalg::Vdense>::null
                                            : value_e),
                          d_value_d_dof_e_flags, d_value_d_dof_e,
                          (d_value_d_time.is_null()
                                 ? b2linalg::Vector<double, b2linalg::Vdense>::null
                                 : d_value_d_time_e));
                }
                if (!value_.is_null()) { value_(index_e) += value_e; }
                if (!d_value_d_dof.is_null()) {
                    d_value_d_dof += b2linalg::StructuredMatrix(
                          d_value_d_dof.size1(), index_e, d_value_d_dof_e[0]);
                }
            }
        }
    }
 
    typedef ObjectTypeComplete<
          NBC_NLR_Coupling, TypedNaturalBoundaryConditionComponent<double>::type_t>
          type_t;
    static type_t type;
 
private:
    using global_elements_t = std::set<TypedElement<double>*>;
    friend class SolutionStaticNonlinearNLRCoupling;
    SolutionStaticNonlinearNLRCoupling solution;
 
    std::string id;
 
    // the local case iterator
    b2000::CaseList::case_iterator local_case_iterator;
 
    // The current local case
    b2000::Case* local_case;
 
    // The fine scale model
    b2000::b2dbv3::Model local_model;
 
    // The coarse scale model
    b2000::Model* global_model;
 
    int subcase_id;
 
    // The interpolation of the coarse scale displacement d2 in the fine scale model
    // N(i, j) * d2[j] is the contribution of the coarse scale dof j in the fine scale dof i
    b2linalg::Matrix<double> C;
    b2linalg::Index C_index;
    bool C_linear;
    std::map<size_t, size_t> C_index_set;
 
    // update the C matrix
    // this is independent of the dof, to call each time the mesh of the local model is updated
    void update_local_model() {
        local_model.init(id);
 
        local_case_iterator = local_model.get_case_list().get_case_iterator();
        local_case = local_case_iterator->next();
        local_model.set_case(*local_case);
 
        solution.set_nlr_coupling(*this);
        dynamic_cast<b2000::b2dbv3::SolutionStaticNonlinear<double>&>(global_model->get_solution())
              .add_sub_solution(solution);
 
        // computation of the global element overlapped by local model
        b2000::Domain& local_domain = local_model.get_domain();
        b2000::Domain& global_domain = global_model->get_domain();
        {
            global_elements.clear();
            std::unique_ptr<b2000::Domain::ElementIterator> local_element_iter =
                  local_domain.get_element_iterator();
            Element* element = nullptr;
            for (element = local_element_iter->next(); element != nullptr;
                 element = local_element_iter->next()) {
                Element* global_element = local_domain.get_parent_element_and_nodes_mapping(
                      global_domain, *element, b2linalg::Matrix<double>::null);
                TypedElement<double>& global_element_typed =
                      dynamic_cast<TypedElement<double>&>(*global_element);
                global_elements.insert(&global_element_typed);
            }
        }
 
        // computation of the kinematic coupling operator from the global dof to the local dof
        // we assume that the coupling operator is linear. Thus this will don't
        // work with structured element (shell or beam).
        {
            C_linear = true;
            // iteration on the local node
            std::unique_ptr<b2000::Domain::NodeIterator> local_node_iter =
                  local_domain.get_node_iterator();
            const std::string field_name = "DISP";
            for (Node* local_node = local_node_iter->next(); local_node != nullptr;
                 local_node = local_node_iter->next()) {
                b2linalg::Vector<double> global_icoor;
                Element* global_element = local_domain.get_parent_element_mapping(
                      global_domain, *local_node, global_icoor);
                b2linalg::Index dof_numbering;
                b2linalg::Matrix<double, b2linalg::Mrectangle> d_value_d_dof;
                TypedElement<double>& e = dynamic_cast<TypedElement<double>&>(*global_element);
#ifdef CHECK_GLOBAL_TO_LOCAL_MAPPING
                {
                    b2linalg::Vector<double> global_coor;
                    e.body_geom(
                          global_icoor, global_coor,
                          b2linalg::Matrix<double, b2linalg::Mrectangle>::null);
                    std::pair<int, const coor_type*> local_coor = local_node->get_coor();
                    for (int i = 0; i != local_coor.first; ++i) {
                        if (norm(local_coor.second[i] - global_coor[i])
                            > CHECK_GLOBAL_TO_LOCAL_MAPPING) {
                            Exception()
                                  << "The local mapping for the local node "
                                  << local_node->get_object_name() << " is not correct." << THROW;
                        }
                    }
                }
#endif
                if (!e.body_field_linear_on_dof(field_name)) { C_linear = false; }
                e.body_field_value(
                      field_name, global_icoor,
                      b2linalg::Matrix<double, b2linalg::Mrectangle_constref>::null, 1,
                      b2linalg::Vector<double, b2linalg::Vdense>::null,
                      b2linalg::Matrix<double, b2linalg::Mrectangle>::null, dof_numbering,
                      d_value_d_dof);
                for (size_t i = 0; i != dof_numbering.size(); ++i) {
                    C_index_set[dof_numbering[i]] = 0;
                }
            }
 
            C_index.resize(C_index_set.size());
            {
                size_t ii = 0;
                for (std::map<size_t, size_t>::iterator i = C_index_set.begin();
                     i != C_index_set.end(); ++i, ++ii) {
                    C_index[ii] = i->first;
                    i->second = ii;
                }
            }
 
            C.resize(local_domain.get_number_of_dof(), C_index.size());
 
            if (C_linear) {
                C.set_zero();
                // iteration on the local node
                get_global_to_local_dof(
                      b2linalg::Vector<double, b2linalg::Vdense_constref>::null,
                      b2linalg::Vector<double, b2linalg::Vdense>::null, C);
            }
        }
 
        // resize the matrix and vector used in the local model
        local_model_residue.init(local_model, *local_case);
        conv_local_dof.clear();
        conv_local_dof.resize(local_model_residue.get_number_of_dof());
        prev_local_dof = conv_local_dof;
    }
 
    void get_global_to_local_dof(
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& dof_global,
          b2linalg::Vector<double, b2linalg::Vdense>& dof_local,
          b2linalg::Matrix<double, b2linalg::Mrectangle_ref> C) {
        b2000::Domain& local_domain = local_model.get_domain();
        b2000::Domain& global_domain = global_model->get_domain();
        const std::string field_name = "DISP";
        if (!dof_local.is_null()) {
            dof_local.set_zero();
            dof_local.resize(local_domain.get_number_of_dof());
        }
        if (!C.is_null()) { C.set_zero(); }
        std::unique_ptr<b2000::Domain::NodeIterator> local_node_iter =
              local_domain.get_node_iterator();
        size_t ii = 0;
        for (Node* local_node = local_node_iter->next(); local_node != nullptr;
             local_node = local_node_iter->next()) {
            b2linalg::Vector<double> global_icoor;
            Element* global_element =
                  local_domain.get_parent_element_mapping(global_domain, *local_node, global_icoor);
            b2linalg::Index dof_numbering;
            b2linalg::Vector<double, b2linalg::Vdense> value;
            b2linalg::Matrix<double, b2linalg::Mrectangle> d_value_d_dof;
            dynamic_cast<TypedElement<double>&>(*global_element)
                  .body_field_value(
                        field_name, global_icoor,
                        (dof_global.is_null()
                               ? b2linalg::Matrix<double, b2linalg::Mrectangle_constref>::null
                               : b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(
                                     dof_global)),
                        1,
                        (dof_local.is_null() ? b2linalg::Vector<double, b2linalg::Vdense>::null
                                             : value),
                        b2linalg::Matrix<double, b2linalg::Mrectangle>::null, dof_numbering,
                        (C.is_null() ? b2linalg::Matrix<double, b2linalg::Mrectangle>::null
                                     : d_value_d_dof));
            if (!dof_local.is_null()) {
                for (size_t i = 0; i != value.size(); ++i) { dof_local[ii + i] = value[i]; }
            }
            if (!C.is_null()) {
#if BCNC_CHECK_DERIVATIVE
                if (!dof_global.is_null()) {
                    const double h = 1e-4;
                    const double inv_12h = 1.0 / (12 * h);
                    b2linalg::Vector<double> dof_h = dof_global;
                    b2linalg::Vector<double> diff;
                    b2linalg::Vector<double> value_tmp;
                    TypedElement<double>& e = dynamic_cast<TypedElement<double>&>(*global_element);
                    bool first_error = true;
                    for (size_t j = 0; j != dof_numbering.size(); ++j) {
                        const size_t jj = dof_numbering[j];
                        dof_h[jj] = dof_global[jj] - 2 * h;
                        e.body_field_value(
                              field_name, global_icoor,
                              b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(dof_h), 1,
                              value_tmp, b2linalg::Matrix<double, b2linalg::Mrectangle>::null,
                              b2linalg::Index::null,
                              b2linalg::Matrix<double, b2linalg::Mrectangle>::null);
                        diff = inv_12h * value_tmp;
 
                        dof_h[jj] = dof_global[jj] + 2 * h;
                        e.body_field_value(
                              field_name, global_icoor,
                              b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(dof_h), 1,
                              value_tmp, b2linalg::Matrix<double, b2linalg::Mrectangle>::null,
                              b2linalg::Index::null,
                              b2linalg::Matrix<double, b2linalg::Mrectangle>::null);
                        diff += -inv_12h * value_tmp;
 
                        dof_h[jj] = dof_global[jj] - h;
                        e.body_field_value(
                              field_name, global_icoor,
                              b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(dof_h), 1,
                              value_tmp, b2linalg::Matrix<double, b2linalg::Mrectangle>::null,
                              b2linalg::Index::null,
                              b2linalg::Matrix<double, b2linalg::Mrectangle>::null);
                        diff += (-8.0 * inv_12h) * value_tmp;
 
                        dof_h[jj] = dof_global[jj] + h;
                        e.body_field_value(
                              field_name, global_icoor,
                              b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(dof_h), 1,
                              value_tmp, b2linalg::Matrix<double, b2linalg::Mrectangle>::null,
                              b2linalg::Index::null,
                              b2linalg::Matrix<double, b2linalg::Mrectangle>::null);
                        diff += (8.0 * inv_12h) * value_tmp;
 
                        dof_h[jj] = dof_global[jj];
 
                        for (size_t i = 0; i != diff.size(); ++i) {
                            if (std::abs(diff[i] - d_value_d_dof(i, j)) > 1e-6) {
                                if (first_error) {
                                    std::cerr << "element " << e.get_object_name()
                                              << ", local node=" << local_node->get_object_name()
                                              << ", coor=" << global_icoor[0] << ", "
                                              << global_icoor[1] << ", " << global_icoor[2]
                                              << std::endl;
                                    first_error = false;
                                }
                                std::cerr << "    d_value_d_dof(" << i << ", " << j
                                          << ")=" << d_value_d_dof(i, j) << ", numeric=" << diff[i]
                                          << std::endl;
                            }
                        }
                    }
                }
#endif
                for (size_t j = 0; j != d_value_d_dof.size2(); ++j) {
                    const size_t jj = C_index_set[dof_numbering[j]];
                    for (size_t i = 0; i != d_value_d_dof.size1(); ++i) {
                        C(ii + i, jj) = d_value_d_dof(i, j);
                    }
                }
            }
            ii += dof_local.is_null() ? d_value_d_dof.size1() : value.size();
        }
    }
 
    // the global dof of the previous iteration (used to compute the delta dof)
    b2linalg::Vector<double> conv_global_dof;
    b2linalg::Vector<double> prev_global_dof;
    b2linalg::Vector<double> conv_local_dof;
    b2linalg::Vector<double> prev_local_dof;
    b2linalg::Matrix<double, b2linalg::Msym_compressed_col_update_inv> local_RtKR;
    b2linalg::Matrix<double> local_RtKC;
 
    global_elements_t global_elements;
 
    StaticResidueFunctionNLRCoupling local_model_residue;
};
 
inline void SolutionStaticNonlinearNLRCoupling::set_nlr_coupling(NBC_NLR_Coupling& nlr_coupling_) {
    nlr_coupling = &nlr_coupling_;
    local_solution = &dynamic_cast<b2000::b2dbv3::SolutionStaticNonlinear<double>&>(
          nlr_coupling->local_model.get_solution());
}
 
inline void SolutionStaticNonlinearNLRCoupling::set_dof(
      double time, const b2linalg::Vector<double, b2linalg::Vdense_constref>& global_dof_all,
      const b2linalg::Vector<double, b2linalg::Vdense_constref>& nbc,
      const b2linalg::Vector<double, b2linalg::Vdense_constref>& residue,
      bool unconverged_subcycle) {
    // compute the local dof without the local reduction
    b2000::Domain& local_domain = nlr_coupling->local_model.get_domain();
    b2linalg::Vector<double, b2linalg::Vdense> g;
    g.resize(local_domain.get_number_of_dof());
    nlr_coupling->conv_local_dof = nlr_coupling->prev_local_dof;
    nlr_coupling->local_model_residue.get_solution(nlr_coupling->conv_local_dof, time, true, g);
 
    // store the local hierarchical dof in the local database
    {
        SetName name = nlr_coupling->local_case->get_string("DOF_SOL", "DISP") + "_H1";
        if (unconverged_subcycle) {
            name.cycle(local_solution->get_cycle_id() + 1);
            name.subcycle(local_solution->get_subcycle_id() + 1);
        } else {
            name.cycle(local_solution->get_cycle_id());
        }
        if (name.case_undef()) { name.case_(nlr_coupling->local_case->get_id()); }
        name.subcase(subcase_id);
        dynamic_cast<DataBaseFieldSet&>(nlr_coupling->local_model)
              .add_field_entry(
                    "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
        RTable rtable;
        rtable.set("STAGE_ID", nlr_coupling->local_case->get_stage_id());
        rtable.set("STAGE", nlr_coupling->local_case->get_stage_number());
        rtable.set("LAMBDA", time);
        dynamic_cast<b2dbv3::Domain&>(nlr_coupling->local_model.get_domain())
              .save_global_dof(
                    name, b2linalg::Vector<double, b2linalg::Vdense_constref>(g), rtable);
    }
 
    // add the global dof to the local hierarchical dof
    {
        if (nlr_coupling->C_linear) {
            b2linalg::Vector<double> global_dof = global_dof_all(nlr_coupling->C_index);
            g += nlr_coupling->C * global_dof;
        } else {
            b2linalg::Vector<double> global_dof_on_local;
            nlr_coupling->get_global_to_local_dof(
                  global_dof_all, global_dof_on_local, nlr_coupling->C);
            g += global_dof_on_local;
        }
    }
 
    // store the local (non hierarchical) dof in the database
    local_solution->set_dof(
          time, g, nlr_coupling->local_model_residue.get_nbc_sol(),
          nlr_coupling->local_model_residue.get_residuum_sol(), unconverged_subcycle);
 
    // store the local (non hierarchical) gradient in the database
    b2000::SolutionStaticNonlinear<double>::gradient_container gc =
          local_solution->get_gradient_container();
    {
        local_domain.set_dof(g);
        b2000::Domain::element_iterator ii = local_domain.get_element_iterator();
        b2linalg::Index index;
        for (Element* e = ii->next(); e != nullptr; e = ii->next()) {
            if (auto te = dynamic_cast<TypedElement<double>*>(e); te != nullptr) {
                te->get_value(
                      nlr_coupling->local_model,
                      false,  // linear
                      true,   // equilibrium_solution
                      time,
                      0,  // delta_time
                      b2linalg::Matrix<double, b2linalg::Mrectangle_constref>(g),
                      gc.get() != nullptr && gc->set_current_element(*e) ? gc.get() : nullptr,
                      nullptr,  // solver_hints
                      index, b2linalg::Vector<double, b2linalg::Vdense>::null,
                      TypedElement<double>::d_value_d_dof_flags_null,
                      TypedElement<double>::d_value_d_dof_null,
                      b2linalg::Vector<double, b2linalg::Vdense>::null);
            }
        }
    }
}
 
inline void SolutionStaticNonlinearNLRCoupling::terminate_stage(bool success) {
    local_solution->terminate_stage(success);
    if (nlr_coupling->local_case->next_stage()) {
        nlr_coupling->local_model.set_case(*nlr_coupling->local_case);
    }
}
 
inline void SolutionStaticNonlinearNLRCoupling::terminate_case(
      bool success, const RTable& attributes) {
    local_solution->terminate_case(success, attributes);
    nlr_coupling->local_case = nlr_coupling->local_case_iterator->next();
    if (nlr_coupling->local_case) { nlr_coupling->local_model.set_case(*nlr_coupling->local_case); }
}
 
}  // namespace b2000::b2dbv3
 
#endif