b2beam_cross_section_solver.H

//------------------------------------------------------------------------
// b2beam_cross_section_solver.H --
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2011-2012,2016 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 B2BEAM_CROSS_SECTION_SOLVER_H_
#define B2BEAM_CROSS_SECTION_SOLVER_H_
 
// #undef DEBUG
 
#include "b2blaslapack.H"
#include "b2ppconfig.h"
#include "elements/smr/b2element_beam_cross_section.H"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model/b2model.H"
#include "model/b2node.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
#include "utils/b2tensor_calculus.H"
#include "utils/b2util.H"
 
namespace b2000::solver {
 
// Reference:
// On the structural behavior and the Saint Venant solution in the exact beam theory:
// Application to laminated composite beams
// Rached El Fatmi and Hatem Zenzri
// Computers and Structures 80 (2002) 1441-1456
 
// #define BCSC_ORIG_BC
// #define BCSC_V1
 
class BeamCrossSectionSolver : public Solver {
public:
    using mrbc_t = TypedModelReductionBoundaryCondition<double>;
    using type_t = ObjectTypeComplete<BeamCrossSectionSolver, Solver::type_t>;
 
    void solve() override {
        while (solve_iteration()) { ; }
    }
 
    bool solve_iteration() override {
        logging::Logger& logger = logging::get_logger("solver.beam_cross_section");
        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_) { return false; }
 
        if (case_->get_stage_size() != 1.0) {
            logging::get_logger("solver.beam_cross_section")
                  << logging::warning << "The stage size of the case " << case_->get_id()
                  << " is equal to " << case_->get_stage_size()
                  << " but the beam cross section solver only handle a case with one stage of size "
                     "equal to one. The specified stage size is ignored."
                  << logging::LOGFLUSH;
        }
        if (case_->get_number_of_stage() != 1) {
            logging::get_logger("solver.beam_cross_section")
                  << logging::warning << "The number of stage of the case" << case_->get_id()
                  << " is " << case_->get_stage_size()
                  << " but the beam cross section solver only handle a case with one stage of size "
                     "equal to one. The extra stages are ignored."
                  << logging::LOGFLUSH;
        }
 
        model_->set_case(*case_);
        Domain& domain = model_->get_domain();
        SolutionBeamCrossSection* solution =
              &dynamic_cast<SolutionBeamCrossSection&>(model_->get_solution());
        if (solution == nullptr) {
            TypeError() << "The solver only accept a solution of type SolutionBeamCrossSection"
                        << THROW;
        }
 
        const size_t number_of_dof = 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<double>())));
        b2linalg::Vector<double, b2linalg::Vdense> r(number_of_dof);
        b2linalg::Matrix<double, b2linalg::Mcompressed_col> trans_R;
        mrbc->get_linear_value(r, trans_R);
 
        {
            // Add constraint to obtain a determined system
            // first find two nodes with different coordinate.
 
            Node* n1 = nullptr;
            bool n1_lookdof1 = false;
            std::pair<int, const double*> coor1;
            {
                Domain::node_iterator i = domain.get_node_iterator();
                for (Node* n = i->next(); n != nullptr; n = i->next()) {
                    coor1 = n->get_coor<double>();
                    if (coor1.first == 3 && trans_R.get_nb_nonzero(coor1.second[0]) == 0) {
                        n1 = n;
                        n1_lookdof1 = true;
                        break;
                    }
                }
            }
            if (n1 == nullptr) {
                Domain::node_iterator i = domain.get_node_iterator();
                for (Node* n = i->next(); n != nullptr; n = i->next()) {
                    coor1 = n->get_coor<double>();
                    if (coor1.first == 3) {
                        n1 = n;
                        break;
                    }
                }
            }
            if (n1 == nullptr) { Exception() << THROW; }
            Node* n2 = nullptr;
            double dist = 0;
            {
                Domain::node_iterator i = domain.get_node_iterator();
                for (Node* n = i->next(); n != nullptr; n = i->next()) {
                    std::pair<int, const double*> coor2 = n->get_coor<double>();
                    if (coor2.first == 3) {
                        double tmp[3];
                        sub_3(coor1.second, coor2.second, tmp);
                        double dtmp = dot_3(tmp, tmp);
                        if (dtmp > dist) {
                            dist = dtmp;
                            n2 = n;
                        }
                    }
                }
            }
            if (n2 == nullptr) { Exception() << THROW; }
 
            b2linalg::Index ind(5);
            n2->get_global_dof_numbering(&ind[2]);
            n1->get_global_dof_numbering(&ind[0]);
            if (n1_lookdof1) { ind.erase(ind.begin()); }
 
            std::set<size_t> col_to_remove;
            for (size_t i = 0; i != ind.size(); ++i) {
                size_t* rowind;
                double* value;
                const size_t nbz = trans_R.get_col_value(ind[i], rowind, value);
                if (nbz == 0) { Exception() << THROW; }
                col_to_remove.insert(rowind, rowind + nbz);
            }
            ind.assign(col_to_remove.begin(), col_to_remove.end());
            trans_R.remove_row(ind);
        }
 
        double neutral_point[2];
        double mass;
        double inertia_point[2];
        double inertia_moment[3];
        b2linalg::Matrix<double> u;
        b2linalg::Matrix<double> lambda;
        b2linalg::Vector<double> lambdad;
        b2linalg::Matrix<double> strain[7];
        b2linalg::Matrix<double> stress[7];
 
        solve_static(
              logger, *case_, model_, domain, trans_R, neutral_point, mass, inertia_point,
              inertia_moment, u, lambda, lambdad, strain, stress);
 
        solution->set_solution(
              neutral_point, mass, inertia_point, inertia_moment, u, strain, stress, lambda,
              lambdad);
 
        RTable attributes;
        solution->terminate_case(true, attributes);
        case_->warn_on_non_used_key();
 
        return true;
    }
 
    static void solve_static(
          logging::Logger& logger, Case& case_, Model* model, Domain& domain,
          b2linalg::Matrix<double, b2linalg::Mcompressed_col>& trans_R, double neutral_point[2],
          double& inertia_mass, double inertia_point[2], double inertia_moment[3],
          b2linalg::Matrix<double>& u, b2linalg::Matrix<double>& lambda,
          b2linalg::Vector<double>& lambdad, b2linalg::Matrix<double> strain[7],
          b2linalg::Matrix<double> stress[7]);
 
    static type_t type;
};
 
}  // namespace b2000::solver
 
// #define TEST_WITOUT_SYMMETRISATION
void b2000::solver::BeamCrossSectionSolver::solve_static(
      logging::Logger& logger, Case& case_, Model* model, Domain& domain,
      b2linalg::Matrix<double, b2linalg::Mcompressed_col>& trans_R, double neutral_point[2],
      double& inertia_mass, double inertia_point[2], double inertia_moment[3],
      b2linalg::Matrix<double>& u, b2linalg::Matrix<double>& lambda,
      b2linalg::Vector<double>& lambdad, b2linalg::Matrix<double> strain[7],
      b2linalg::Matrix<double> stress[7]) {
    logger << logging::info << "Start the beam cross section solver" << logging::LOGFLUSH;
 
    const size_t number_of_dof = domain.get_number_of_dof();
 
    const bool log_el_flag = logger.is_enabled_for(logging::data);
    logger << logging::info << "Computation of centre, first and second moment of area"
           << logging::LOGFLUSH;
 
    neutral_point[0] = neutral_point[1] = 0;
    inertia_mass = 0;
    inertia_point[0] = inertia_point[1] = 0;
    {
        double area = 0;
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_0(
                  model, area, neutral_point, inertia_mass, inertia_point);
        }
        neutral_point[0] /= area;
        neutral_point[1] /= area;
        if (inertia_mass > 0) {
            inertia_point[0] = inertia_point[0] / inertia_mass - neutral_point[0];
            inertia_point[1] = inertia_point[1] / inertia_mass - neutral_point[1];
        }
    }
 
    double inertia_axis[2][2];
    {
        double inertia[3] = {};
        std::fill_n(inertia_moment, 3, 0);
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_01(
                  model, neutral_point, inertia, inertia_moment);
        }
        {
            double value[2];
            double m[2][2] = {{inertia[0], inertia[2]}, {inertia[2], inertia[1]}};
            eigenvector_2(m, value, inertia_axis);
            normalise_2(inertia_axis[0]);
            normalise_2(inertia_axis[1]);
            if (determinant_2_2(inertia_axis) < 0) {
                for (int j = 0; j != 2; ++j) { inertia_axis[1][j] *= -1; }
            }
        }
    }
 
    logger << logging::info << "Element matrix assembly" << logging::LOGFLUSH;
 
    // we add the diagonal one for additional dof
    for (int i = 0; i != 4; ++i) { trans_R.push_back(trans_R.size1(), trans_R.size2(), 1); }
 
    size_t K_augm_size = trans_R.size1();
    b2linalg::Matrix<double, b2linalg::Msym_compressed_col_update_inv> K_augm(
          K_augm_size, domain.get_dof_connectivity_type(), case_);
    domain.set_dof(b2linalg::Vector<double, b2linalg::Vdense_constref>::null);
    {
        Domain::element_iterator i = domain.get_element_iterator();
        b2linalg::Index index;
        b2linalg::Matrix<double, b2linalg::Mpacked> Ke;
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_1(
                  model, neutral_point, inertia_axis, number_of_dof, index, Ke);
            K_augm +=
                  trans_R * StructuredMatrix(number_of_dof + 4, index, Ke) * transposed(trans_R);
            if (log_el_flag) {
                b2linalg::Matrix<double> tmp(Ke);
                logger << logging::data << "elementary second variation "
                       << logging::DBName("ELSV", 0, e->get_object_name()) << tmp
                       << " of element id " << e->get_object_name() << logging::LOGFLUSH;
            }
        }
    }
 
    logger << logging::info << "First system resolution" << logging::LOGFLUSH;
 
    b2linalg::Matrix<double> v(number_of_dof + 4, 6);
    lambda.resize(6, 6);
    {
        b2linalg::Matrix<double> v_augm(K_augm.size1(), 4);
        for (int i = 0; i != 4; ++i) { v_augm(K_augm.size1() - 4 + i, i) = 1; }
        v_augm = b2linalg::inverse(K_augm) * v_augm;
        v(b2linalg::Interval(2, 6)) = transposed(trans_R) * v_augm;
        v[0] = v[2];
    }
    for (int j = 0; j != 4; ++j) {
        const int jj = j == 0 ? 0 : j + 2;
        for (int i = 0; i != 4; ++i) {
            const int ii = i == 0 ? 0 : i + 2;
#ifdef TEST_WITOUT_SYMMETRISATION
            lambda(ii, jj) = v(number_of_dof + i, jj);
#else
            lambda(ii, jj) = 0.5 * (v(number_of_dof + i, jj) + v(number_of_dof + j, ii));
#endif
        }
    }
 
    logger << logging::info << "Factorisation and first computation of beta" << logging::LOGFLUSH;
 
    b2float128 beta[6][6] = {};
    {
        const bool beta_cd[6][6] = {{1, 0, 0, 1, 1, 1}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0},
                                    {1, 0, 0, 1, 1, 1}, {1, 0, 0, 1, 1, 1}, {1, 0, 0, 1, 1, 1}};
        const int beta_c[] = {0, 3, 4, 5};
        const int beta_d[] = {0, 3, 4, 5};
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_2(
                  model, neutral_point, inertia_axis, v, lambda, beta_cd, 4, beta_c, 4, beta_d,
                  beta);
        }
    }
    for (int j = 0; j != 4; ++j) {
        const int jj = j == 0 ? 0 : j + 2;
#ifdef TEST_WITOUT_SYMMETRISATION
        lambda(1, jj) = beta[jj][5];
        lambda(2, jj) = -beta[jj][4];
#else
        lambda(1, jj) = lambda(jj, 1) = beta[jj][5];
        lambda(2, jj) = lambda(jj, 2) = -beta[jj][4];
#endif
    }
 
    logger << logging::info << "Second system resolution" << logging::LOGFLUSH;
    {
        b2linalg::Matrix<double> u(number_of_dof + 4, 2);
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_3(
                  model, neutral_point, inertia_axis, number_of_dof, v, lambda, beta, u);
        }
 
        b2linalg::Matrix<double> u_augm;
        u_augm = trans_R * u;
        u_augm = b2linalg::inverse(K_augm) * u_augm;
        v(b2linalg::Interval(1, 3)) = transposed(trans_R) * u_augm;
    }
    for (int j = 1; j != 3; ++j) {
        for (int i = 0; i != 4; ++i) {
            const int ii = i == 0 ? 0 : i + 2;
#ifdef TEST_WITOUT_SYMMETRISATION
            lambda(ii, j) = v(number_of_dof + i, j);
#else
            lambda(ii, j) = 0.5 * (lambda(ii, j) + v(number_of_dof + i, j));
            lambda(j, ii) = 0.5 * (lambda(j, ii) + v(number_of_dof + i, j));
#endif
        }
    }
 
    // computation of u from v
    u.resize(number_of_dof, 7);
    for (int i = 0; i != 6; ++i) {
        const double tmp1[3] = {double(beta[i][0]), double(beta[i][1]), double(beta[i][2])};
        double tmp2[3] = {double(beta[i][3]), double(beta[i][4]), double(beta[i][5])};
 
        Domain::node_iterator iter = domain.get_node_iterator();
        for (Node* n = iter->next(); n != nullptr; n = iter->next()) {
            double tmp3[3];
            const double* ptr = n->get_coor<double>().second;
            const double tmp4[3] = {0, ptr[1] - neutral_point[0], ptr[2] - neutral_point[1]};
            outer_product_3(tmp2, tmp4, tmp3);
            const size_t pos = n->get_global_dof_numbering().first;
            for (int ii = 0; ii != 3; ++ii) {
                u(pos + ii, i) = v(pos + ii, i) + tmp1[ii] + tmp3[ii];
            }
        }
    }
 
    logger << logging::info << "Second computation of beta" << logging::LOGFLUSH;
    {
        const bool beta_cd[6][6] = {{0, 1, 1, 0, 0, 0}, {1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1},
                                    {0, 1, 1, 0, 0, 0}, {0, 1, 1, 0, 0, 0}, {0, 1, 1, 0, 0, 0}};
        const int beta_c[] = {0, 1, 2, 3, 4, 5};
        const int beta_d[] = {0, 1, 2, 3, 4, 5};
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_2(
                  model, neutral_point, inertia_axis, v, lambda, beta_cd, 6, beta_c, 6, beta_d,
                  beta);
        }
    }
    lambda(1, 1) = beta[5][1] + beta[1][5];
    lambda(2, 1) = beta[5][2] - beta[1][4];
    lambda(1, 2) = -beta[4][1] + beta[2][5];
    lambda(2, 2) = -beta[4][2] - beta[2][4];
#ifndef TEST_WITOUT_SYMMETRISATION
    lambda(2, 1) = 0.5 * (lambda(2, 1) + lambda(1, 2));
    lambda(1, 2) = lambda(2, 1);
#endif
 
    logger << logging::info << "Strain and stress computation" << logging::LOGFLUSH;
 
    const size_t number_of_node = domain.get_number_of_nodes();
    {
        for (int i = 0; i != 7; ++i) {
            strain[i].resize(6, number_of_node);
            stress[i].resize(3, number_of_node);
        }
    }
 
    std::vector<size_t> node_weight(number_of_dof);
    std::vector<std::pair<const Node*, const Node*> > slave_nodes;
    domain.get_slave_node(slave_nodes);
    {
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_4(
                  model, neutral_point, inertia_axis, v, lambda, beta, node_weight, strain, stress);
        }
 
        for (size_t ptr = 0; ptr != slave_nodes.size(); ++ptr) {
            const size_t sn = slave_nodes[ptr].first->get_id();
            const size_t mn = slave_nodes[ptr].second->get_id();
            node_weight[mn] += node_weight[sn];
            for (int i = 0; i != 6; ++i) {
                for (int j = 0; j != 6; ++j) { strain[i](j, mn) += strain[i](j, sn); }
                for (int j = 0; j != 3; ++j) { stress[i](j, mn) += stress[i](j, sn); }
            }
        }
 
        for (size_t k = 0; k != number_of_node; ++k) {
            const double s = 1.0 / node_weight[k];
            for (int i = 0; i != 6; ++i) {
                for (int j = 0; j != 6; ++j) { strain[i](j, k) *= s; }
                for (int j = 0; j != 3; ++j) { stress[i](j, k) *= s; }
            }
        }
 
        for (size_t ptr = 0; ptr != slave_nodes.size(); ++ptr) {
            const size_t sn = slave_nodes[ptr].first->get_id();
            const size_t mn = slave_nodes[ptr].second->get_id();
            for (int i = 0; i != 6; ++i) {
                for (int j = 0; j != 6; ++j) { strain[i](j, sn) = strain[i](j, mn); }
                for (int j = 0; j != 3; ++j) { stress[i](j, sn) = stress[i](j, mn); }
            }
        }
    }
 
    logger << logging::info << "Computation of f due to the initial stress and bc"
           << logging::LOGFLUSH;
    double f[6];
    std::fill_n(f, 6, 0);
    {
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_5(
                  model, neutral_point, inertia_axis, f);
        }
    }
 
    b2linalg::Vector<double> vd(trans_R.size2());
#ifdef BCSC_V1
    logger << logging::info << "Third system resolution due to initial stress and bc"
           << logging::LOGFLUSH;
    {
        b2linalg::Vector<double> value(number_of_dof + 4);
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != 0; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_6(
                  model, neutral_point, inertia_axis, number_of_dof, f, u, stress, value);
        }
#ifdef BCSC_ORIG_BC
        b2linalg::Vector<double> u_augm;
        u_augm = trans_R * value;
        u_augm = b2linalg::inverse(K_augm) * u_augm;
        vd = transposed(trans_R) * u_augm;
 
        lambdad.resize(6);
        lambdad[0] = vd[number_of_dof];
        for (int i = 1; i != 4; ++i) { lambdad[i + 2] = vd[number_of_dof + i]; }
#else
        b2linalg::Matrix<double> u_augm(trans_R.size1(), 5);
        u_augm[4] = trans_R * value;
        for (int i = 0; i != 4; ++i) { u_augm(trans_R.size1() - 4 + i, i) = 1; }
        u_augm = b2linalg::inverse(K_augm) * u_augm;
 
        b2linalg::Vector<double> lm;
        b2linalg::Interval int1(K_augm.size1() - 4, K_augm.size1());
        // lm = b2linalg::inverse(u_augm(int1, b2linalg::Interval(0, 4))) * u_augm[4](int1);
        {
            int ipiv[4];
            int info;
            lm = u_augm[4](int1);
            lapack::gesv(
                  4, 1, &u_augm(K_augm.size1() - 4, 0), u_augm.size1(), ipiv, &lm[0], 4, info);
            if (info != 0) { Exception() << THROW; }
        }
 
        u_augm[4] = trans_R * value;
        u_augm[4](int1) -= lm;
        {
            b2linalg::Vector<double> tmp;
            tmp = b2linalg::inverse(K_augm) * u_augm[4];
            vd = transposed(trans_R) * tmp;
        }
 
        lambdad.resize(6);
        lambdad[0] = lm[0];
        for (int i = 1; i != 4; ++i) { lambdad[i + 2] = lm[i]; }
#endif
    }
#else
    lambdad.resize(6);
    for (int i = 0; i != 6; ++i) { lambdad[i] = f[i]; }
#endif
 
    logger << logging::info << "Strain and stress computation due to initial stress and bc"
           << logging::LOGFLUSH;
    {
        b2linalg::Vector<double> value;
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_7(
                  model, neutral_point, inertia_axis, f, u, vd,
                  vd(b2linalg::Interval(number_of_dof, number_of_dof + 4)), strain[6], stress[6]);
        }
 
        for (size_t ptr = 0; ptr != slave_nodes.size(); ++ptr) {
            const size_t sn = slave_nodes[ptr].first->get_id();
            const size_t mn = slave_nodes[ptr].second->get_id();
            for (int j = 0; j != 6; ++j) { strain[6](j, mn) += strain[6](j, sn); }
            for (int j = 0; j != 3; ++j) { stress[6](j, mn) += stress[6](j, sn); }
        }
 
        for (size_t k = 0; k != number_of_node; ++k) {
            const double s = 1.0 / node_weight[k];
            for (int j = 0; j != 6; ++j) { strain[6](j, k) *= s; }
            for (int j = 0; j != 3; ++j) { stress[6](j, k) *= s; }
        }
 
        for (size_t ptr = 0; ptr != slave_nodes.size(); ++ptr) {
            const size_t sn = slave_nodes[ptr].first->get_id();
            const size_t mn = slave_nodes[ptr].second->get_id();
            for (int i = 0; i != 6; ++i) {
                for (int j = 0; j != 6; ++j) { strain[6](j, sn) = strain[6](j, mn); }
                for (int j = 0; j != 3; ++j) { stress[6](j, sn) = stress[6](j, mn); }
            }
        }
    }
 
#ifdef BCSC_V1
    logger << logging::info << "beta computation due to initial stress and bc" << logging::LOGFLUSH;
    {
        double betad[6];
        std::fill_n(betad, 6, 0);
        Domain::element_iterator i = domain.get_element_iterator();
        for (Element* e = i->next(); e != 0; e = i->next()) {
            static_cast<BeamCrossSectionBaseElement&>(*e).add_value_8(
                  model, neutral_point, inertia_axis, u, vd, stress, stress[6], betad);
        }
        lambdad[1] = betad[5];
        lambdad[2] = -betad[4];
 
        // computation of ud from vd
        const double tmp1[3] = {betad[0], betad[1], betad[2]};
        double tmp2[3] = {betad[3], betad[4], betad[5]};
 
        Domain::node_iterator iter = domain.get_node_iterator();
        for (Node* n = iter->next(); n != 0; n = iter->next()) {
            double tmp3[3];
            const double* ptr = n->get_coor<double>().second;
            const double tmp4[3] = {0, ptr[1] - neutral_point[0], ptr[2] - neutral_point[1]};
            outer_product_3(tmp2, tmp4, tmp3);
            const size_t pos = n->get_global_dof_numbering().first;
            for (int ii = 0; ii != 3; ++ii) { u(pos + ii, 6) = vd[pos + ii] + tmp1[ii] + tmp3[ii]; }
        }
    }
#endif
 
    // update u in function of the inertia axis
    {
        for (int i = 0; i != 7; ++i) {
            Domain::node_iterator iter = domain.get_node_iterator();
            for (Node* n = iter->next(); n != nullptr; n = iter->next()) {
                const size_t pos = n->get_global_dof_numbering().first;
                const double tmp[2] = {u(pos + 1, i), u(pos + 2, i)};
                u(pos + 1, i) = inertia_axis[0][0] * tmp[0] + inertia_axis[1][0] * tmp[1];
                u(pos + 2, i) = inertia_axis[0][1] * tmp[0] + inertia_axis[1][1] * tmp[1];
            }
        }
    }
 
    // update the strain and stress in function of the inertia axis
    {
        for (size_t k = 0; k != number_of_node; ++k) {
            for (int i = 0; i != 7; ++i) {
                const double tmp[2] = {stress[i](1, k), stress[i](2, k)};
                stress[i](1, k) = inertia_axis[0][0] * tmp[0] + inertia_axis[1][0] * tmp[1];
                stress[i](2, k) = inertia_axis[0][1] * tmp[0] + inertia_axis[1][1] * tmp[1];
            }
        }
    }
 
    // update lambda in function of the inertia axis
    {
        b2linalg::Matrix<double> trans(6, 6);
        trans(0, 0) = trans(3, 3) = 1;
        for (int j = 0; j != 2; ++j) {
            for (int i = 0; i != 2; ++i) {
                trans(i + 1, j + 1) = trans(i + 4, j + 4) = inertia_axis[j][i];
            }
        }
        /*
        trans(4, 0) = -neutral_point[1];
        trans(5, 0) = neutral_point[0];
        trans(3, 1) = neutral_point[1] * inertia_axis[0][0] - neutral_point[0] * inertia_axis[0][1];
        trans(3, 2) = neutral_point[1] * inertia_axis[1][0] - neutral_point[0] * inertia_axis[1][1];
        */
 
        {
            b2linalg::Matrix<double> tmp;
            tmp = trans * lambda;
            lambda = tmp * transposed(trans);
        }
#ifdef BCSC_ORIG_BC
        {
            b2linalg::Vector<double> tmp;
            tmp = trans * lambdad;
            lambdad = tmp * transposed(trans);
        }
#else
        {
            b2linalg::Vector<double> tmp;
            tmp = transposed(trans) * lambdad;
            lambdad = tmp * trans;
        }
#endif
    }
 
#ifdef DEBUG
    std::cout << "beta=" << std::endl;
    for (int i = 0; i != 6; ++i) {
        for (int j = 0; j != 6; ++j) { std::cout << double(beta[j][i]) << " "; }
        std::cout << std::endl;
    }
    double sum[4] = {};
    for (int i = 0; i != 3; ++i) {
        sum[0] += beta[i][i];
        sum[1] += beta[i + 3][i + 3];
        sum[2] += beta[i][i + 3];
        sum[3] += beta[i + 3][i];
    }
    std::cout << "localisation prop: " << sum[0] << " " << sum[1] << " " << sum[2] << " " << sum[3]
              << std::endl;
    std::cout << "lambda=\n" << lambda << std::endl;
#endif
 
    b2linalg::Matrix<double, b2linalg::Mrectangle_increment_ref>(lambda).in_place_invert();
#ifdef BCSC_ORIG_BC
    {
        b2linalg::Vector<double> tmp;
        tmp = lambda * lambdad;
        lambdad = -tmp;
    }
#endif
 
#ifdef DEBUG
    std::cout << lambda << std::endl;
 
    {
        b2linalg::Matrix<double> lambda_svd(lambda);
        double work[48];
        Vector<double> eigenvalue(6);
        int info;
        lapack::syev('V', 'L', 6, &lambda_svd(0, 0), 6, &eigenvalue[0], work, 48, info);
        assert(info == 0);
        std::cout << "eigenvalue=\n" << eigenvalue << std::endl;
        std::cout << "lambda_svd=\n" << lambda_svd << std::endl;
    }
#endif
}
 
#endif  // B2BEAM_CROSS_SECTION_SOLVER_H_