21#ifndef B2ELEMENT_SHELL_MITC_H_ 
   22#define B2ELEMENT_SHELL_MITC_H_ 
   28#include "b2element_continuum_integrate.H" 
   29#include "b2element_continuum_util.H" 
   30#include "elements/properties/b2solid_mechanics.H" 
   31#include "model/b2element.H" 
   32#include "model_imp/b2hnode.H" 
   34#include "utils/b2linear_algebra.H" 
   35#include "utils/b2rtable.H" 
   39#define TT_INTEGRATION_AT_GAUSS_POINT 3 
   44struct MITCTyingPointDef {
 
   45    static const int nb_tying_point = 0;
 
   53    static const TyingPoint tying_point[nb_tying_point];
 
   55    static const int nb_err_comp = 0;
 
   56    static const int nb_ess_comp = 0;
 
   57    static const int nb_ers_comp = 0;
 
   58    static const int nb_ert_comp = 0;
 
   59    static const int nb_est_comp = 0;
 
   67    struct TyingPointInterp {
 
   68        Interp err[nb_err_comp];
 
   69        Interp ess[nb_ess_comp];
 
   70        Interp ers[nb_ers_comp];
 
   71        Interp ert[nb_ert_comp];
 
   72        Interp est[nb_est_comp];
 
   75    static void get_interpolation(
double r, 
double s, TyingPointInterp& interp) {}
 
   79      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
   80      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
   81      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
   84    using dof2_node_type = node::HNode<
 
   85          coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
 
   86          coordof::Dof<coordof::DTrans<coordof::RX, coordof::RY>>>;
 
   88    using dof5_node_type = node::HNode<
 
   90                coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
 
   91                coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
 
   93                coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
 
   94                coordof::DTrans<coordof::RX, coordof::RY>>>;
 
   96    using dof6_node_type = node::HNode<
 
   98                coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
 
   99                coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
 
  101                coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
 
  102                coordof::DTrans<coordof::RX, coordof::RY, coordof::RZ>>>;
 
  104    using type_t = ObjectTypeComplete<ShellMITC, typename TypedElement<TP>::type_t>;
 
  107        : incompatible_dof_index(0),
 
  109          non_structural_mass(0),
 
  110          line_load_as_moment(false),
 
  111          properties_list(0) {}
 
  114        if (properties_list) {
 
  115            if (property->linear() != SolidMechanics::yes) {
 
  116                const bool shell_interface = 
property->has_shell_stress_interface();
 
  117                int number_of_layer = shell_interface ? 1 : 
property->number_of_layer();
 
  118                for (
int l = 0; l != number_of_layer; ++l) {
 
  119                    for (
int g = 0; g != nb_gauss; ++g) {
 
  120                        for (
int i = 0; i != (shell_interface ? 1 : nb_thickness_gauss); ++i) {
 
  121                            delete properties_list[l][g][i];
 
  126            delete[] properties_list;
 
  131    int get_number_of_dof()
 const override {
 
  132        if (incompatible_mode) {
 
  139    size_t set_global_dof_numbering(
size_t index)
 override {
 
  140        incompatible_dof_index = index;
 
  141        if (incompatible_mode) {
 
  148    std::pair<size_t, size_t> get_global_dof_numbering()
 const override {
 
  149        if (incompatible_mode) {
 
  150            return std::pair<size_t, size_t>(incompatible_dof_index, incompatible_dof_index + 4);
 
  152            return std::pair<size_t, size_t>(0, 0);
 
  156    void set_nodes(std::pair<int, Node* const*> nodes_)
 override {
 
  157        if (nodes_.first != nb_nodes) {
 
  158            Exception() << 
"This element has " << nb_nodes << 
" nodes." << 
THROW;
 
  160        for (
int i = 0; i != nb_nodes_ip; ++i) {
 
  161            nodes[i] = nodes_.second[i];
 
  162            if (!dof5_node_type::is_dof_subset_s(nodes[i])) {
 
  164                            << 
" cannot accept node " << nodes[i]->get_object_name()
 
  165                            << 
" because it does not have the dofs UX,UY,UZ,RX,RY." << 
THROW;
 
  167            nodes_type_6_dof[i] = dof6_node_type::is_dof_subset_s(nodes[i]);
 
  169        for (
int i = nb_nodes_ip; i != nb_nodes; ++i) {
 
  170            if (!dof2_node_type::is_dof_subset_s(nodes[i])) {
 
  172                            << 
" cannot accept node " << nodes[i]->get_object_name()
 
  173                            << 
" because it does not have the dofs RX,RY." << 
THROW;
 
  175            nodes_type_6_dof[i] = 
false;
 
  179    std::pair<int, Node* const*> get_nodes()
 const override {
 
  180        return std::pair<int, Node* const*>(nb_nodes, nodes);
 
  183    void set_property(ElementProperty* property_)
 override {
 
  184        property = 
dynamic_cast<SolidMechanics25D<TP>*
>(property_);
 
  186            Exception() << 
"Bad or missing element property type for shell " 
  188                        << this->
get_object_name() << 
" (forgot MID or PID element attribute?) " 
  193    const ElementProperty* get_property()
 const override { 
return property; }
 
  195    void set_additional_properties(
const RTable& rtable)
 override {
 
  196        non_structural_mass = rtable.get<TP>(
"NON_STRUCTURAL_MASS", TP(0));
 
  197        if (rtable.has_key(
"LINE_LOAD")) {
 
  198            const std::string s = rtable.get_string(
"LINE_LOAD");
 
  199            if (s == 
"MOMENT") { line_load_as_moment = 
true; }
 
  203    void add_dof_field(Element::ListDofField& ldfield)
 override {
 
  204        if (incompatible_mode) {
 
  207                  std::pair<size_t, size_t>(incompatible_dof_index, incompatible_dof_index + 4));
 
  211    void init(Model& model)
 override {
 
  212        if (property->linear() == SolidMechanics::yes) {
 
  215        } 
else if (property->path_dependent()) {
 
  216            const bool shell_interface = 
property->has_shell_stress_interface();
 
  217            int number_of_layer = shell_interface ? 1 : 
property->number_of_layer();
 
  218            if (properties_list) {
 
  219                for (
int l = 0; l != number_of_layer; ++l) {
 
  220                    for (
int g = 0; g != nb_gauss; ++g) {
 
  221                        for (
int i = 0; i != (shell_interface ? 1 : nb_thickness_gauss); ++i) {
 
  222                            delete properties_list[l][g][i];
 
  226                delete[] properties_list;
 
  229                  new SolidMechanics25D<TP>*[number_of_layer][nb_gauss][nb_thickness_gauss];
 
  230            for (
int l = 0; l != number_of_layer; ++l) {
 
  231                for (
int g = 0; g != nb_gauss; ++g) {
 
  232                    for (
int i = 0; i != (shell_interface ? 1 : nb_thickness_gauss); ++i) {
 
  233                        properties_list[l][g][i] =
 
  234                              dynamic_cast<SolidMechanics25D<TP>*
>(
property->copy(l));
 
  247            double R[nb_nodes][6] = {};
 
  248            for (
int k = 0; k != nb_nodes; ++k) {
 
  249                if (nodes_type_6_dof[k]) {
 
  250                    dof6_node_type::get_coor_s(R[k], nodes[k]);
 
  252                    dof5_node_type::get_coor_s(R[k], nodes[k]);
 
  256            const double eps = 1e-12;
 
  257            bool inverted = 
false;
 
  263                b2linalg::Matrix<double> dN(2, nb_nodes);
 
  264                SHAPE::eval_h_derived(xi, dN);
 
  268                for (
int k = 0; k != nb_nodes; ++k) {
 
  269                    for (
int j = 0; j != 3; ++j) {
 
  270                        Ar[j] += dN[k][0] * R[k][j];
 
  271                        As[j] += dN[k][1] * R[k][j];
 
  275                if (
normalise_3(n_center) < eps) { inverted = 
true; }
 
  278            INTEGRATION integration;
 
  279            integration.set_num(nb_gauss);
 
  280            for (
int l = 0; !inverted && l != nb_gauss; ++l) {
 
  282                integration.get_point(l, ip_id, gauss_point[l].weight);
 
  284                b2linalg::Matrix<double> dN(2, nb_nodes);
 
  285                SHAPE::eval_h_derived(ip_id.data(), dN);
 
  291                for (
int k = 0; k != nb_nodes; ++k) {
 
  292                    for (
int j = 0; j != 3; ++j) {
 
  293                        Ar[j] += dN[k][0] * R[k][j];
 
  294                        As[j] += dN[k][1] * R[k][j];
 
  305                if (
dot_3(n, n_center) < eps) {
 
  312                            << 
" is inverted." << 
THROW;
 
  317    void get_dof_numbering(b2linalg::Index& dof_numbering)
 override {
 
  318        dof_numbering.resize(get_number_of_all_dof());
 
  319        size_t* ptr = &dof_numbering[0];
 
  320        for (
int k = 0; k != nb_nodes_ip; ++k) {
 
  321            if (nodes_type_6_dof[k]) {
 
  322                ptr = dof6_node_type::get_global_dof_numbering_s(ptr, nodes[k]);
 
  324                ptr = dof5_node_type::get_global_dof_numbering_s(ptr, nodes[k]);
 
  327        for (
int k = nb_nodes_ip; k != nb_nodes; ++k) {
 
  328            ptr = dof2_node_type::get_global_dof_numbering_s(ptr, nodes[k]);
 
  330        if (incompatible_mode) {
 
  331            for (
int i = 0; i != 4; ++i) { ptr[i] = incompatible_dof_index + i; }
 
  335    const std::vector<Element::VariableInfo> get_value_info()
 const override {
 
  340          Model& model, 
const bool linear, 
const EquilibriumSolution equilibrium_solution,
 
  341          const double time, 
const double delta_time,
 
  342          const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof,
 
  343          GradientContainer* gradient_container, SolverHints* solver_hints,
 
  344          b2linalg::Index& dof_numbering, b2linalg::Vector<TP, b2linalg::Vdense>& value,
 
  345          const std::vector<bool>& d_value_d_dof_flags,
 
  346          std::vector<b2linalg::Matrix<TP, b2linalg::Mpacked>>& d_value_d_dof,
 
  347          b2linalg::Vector<TP, b2linalg::Vdense>& d_value_d_time) 
override;
 
  349    int edge_field_order(
const int edge_id, 
const std::string& field_name)
 override {
 
  350        return SHAPE::edge_shape_type_0::order;
 
  353    bool edge_field_linear_on_dof(
const int edge_id, 
const std::string& field_name)
 override {
 
  357    int edge_field_polynomial_sub_edge(
 
  358          const int edge_id, 
const std::string& field_name,
 
  359          b2linalg::Vector<double, b2linalg::Vdense>& sub_nodes)
 override {
 
  361        sub_nodes.reserve(2);
 
  362        sub_nodes.push_back(-1);
 
  363        sub_nodes.push_back(1);
 
  368          const int edge_id, 
const double internal_coor, b2linalg::Vector<TP>& geom,
 
  369          b2linalg::Vector<TP>& d_geom_d_icoor) {
 
  370        typedef typename SHAPE::edge_shape_type_0 edge_shape;
 
  371        const int* edge_node = SHAPE::edge_node[edge_id - 1];
 
  372        TP node_coor[edge_shape::num_nodes][3] = {};
 
  373        for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  374            if (nodes_type_6_dof[edge_node[k]]) {
 
  375                dof6_node_type::get_coor_s(node_coor[k], nodes[edge_node[k]]);
 
  377                dof5_node_type::get_coor_s(node_coor[k], nodes[edge_node[k]]);
 
  380        if (!geom.is_null()) {
 
  381            b2linalg::Vector<double> shape(edge_shape::num_nodes);
 
  382            edge_shape::eval_h(&internal_coor, shape);
 
  384            for (
size_t i = 0; i != 3; ++i) {
 
  386                for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  387                    tmp += shape[k] * node_coor[k][i];
 
  392        if (!d_geom_d_icoor.is_null()) {
 
  393            b2linalg::Matrix<double> d_shape(1, edge_shape::num_nodes);
 
  394            edge_shape::eval_h_derived(&internal_coor, d_shape);
 
  395            d_geom_d_icoor.resize(3);
 
  396            for (
size_t i = 0; i != 3; ++i) {
 
  398                for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  399                    tmp += d_shape(0, k) * node_coor[k][i];
 
  401                d_geom_d_icoor[i] = tmp;
 
  406    void edge_field_value(
 
  407          const int edge_id, 
const std::string& field_name, 
const double internal_coor,
 
  408          const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof, 
const double time,
 
  409          b2linalg::Vector<TP, b2linalg::Vdense>& value,
 
  410          b2linalg::Vector<TP, b2linalg::Vdense>& d_value_d_icoor, b2linalg::Index& dof_numbering,
 
  411          b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_value_d_dof,
 
  412          b2linalg::Index& d_value_d_dof_dep) {
 
  413        typedef typename SHAPE::edge_shape_type_0 edge_shape;
 
  414        const int* edge_node = SHAPE::edge_node[edge_id - 1];
 
  415        size_t dn[3 * edge_shape::num_nodes + 3];
 
  417            size_t* ptr = &dn[0];
 
  418            for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  419                if (nodes_type_6_dof[edge_node[k]]) {
 
  420                    ptr = dof6_node_type::get_global_dof_numbering_s(ptr, nodes[edge_node[k]]) - 3;
 
  422                    ptr = dof5_node_type::get_global_dof_numbering_s(ptr, nodes[edge_node[k]]) - 2;
 
  426        if (!value.is_null()) {
 
  427            b2linalg::Vector<double> shape(edge_shape::num_nodes);
 
  428            edge_shape::eval_h(&internal_coor, shape);
 
  430            for (
size_t i = 0; i != 3; ++i) {
 
  432                for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  433                    tmp += shape[k] * dof(dn[k * 3 + i], 0);
 
  438        if (!d_value_d_icoor.is_null()) {
 
  439            b2linalg::Matrix<double> d_shape(1, edge_shape::num_nodes);
 
  440            edge_shape::eval_h_derived(&internal_coor, d_shape);
 
  441            d_value_d_icoor.resize(3);
 
  442            for (
size_t i = 0; i != 3; ++i) {
 
  444                for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  445                    tmp += d_shape(0, k) * dof(dn[k * 3 + i], 0);
 
  447                d_value_d_icoor[i] = tmp;
 
  450        if (!d_value_d_dof.is_null()) {
 
  451            dof_numbering.resize(edge_shape::num_nodes * 3);
 
  452            std::copy(dn, dn + edge_shape::num_nodes * 3, dof_numbering.begin());
 
  453            b2linalg::Vector<double> shape(edge_shape::num_nodes);
 
  454            edge_shape::eval_h(&internal_coor, shape);
 
  455            d_value_d_dof.resize(3, edge_shape::num_nodes * 3);
 
  456            d_value_d_dof.set_zero();
 
  457            for (
int k = 0; k != edge_shape::num_nodes; ++k) {
 
  458                for (
size_t i = 0; i != 3; ++i) { d_value_d_dof(i, k * 3 + i) = shape[k]; }
 
  463    int face_field_order(
const int face_id, 
const std::string& field_name) {
 
  467            return SHAPE::edge_shape_type_0::order + 1;
 
  471    bool face_field_linear_on_dof(
const int face_id, 
const std::string& field_name)
 override {
 
  477    int face_field_polynomial_sub_face(
 
  478          const int face_id, 
const std::string& field_name,
 
  479          b2linalg::Matrix<double, b2linalg::Mrectangle>& sub_nodes,
 
  480          std::vector<Element::Triangle>& sub_faces)
 override {
 
  482        sub_nodes.resize(2, 4);
 
  483        sub_nodes(0, 0) = -1;
 
  484        sub_nodes(1, 0) = -1;
 
  486        sub_nodes(1, 1) = -1;
 
  489        sub_nodes(0, 3) = -1;
 
  492        sub_faces.reserve(2);
 
  493        sub_faces.push_back(Element::Triangle(0, 1, 2));
 
  494        sub_faces.push_back(Element::Triangle(2, 3, 0));
 
  495        return face_field_order(face_id, field_name);
 
  500          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
 
  501          b2linalg::Vector<TP>& geom, b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_geom_d_icoor) {
 
  504        typedef typename SHAPE::edge_shape_type_0 edge_shape;
 
  505        TP R[nb_nodes][6] = {};
 
  506        TP D[nb_nodes][3] = {};
 
  508        for (
int k = 0; k != nb_nodes; ++k) {
 
  509            if (nodes_type_6_dof[k]) {
 
  510                dof6_node_type::get_coor_s(R[k], nodes[k]);
 
  512                dof5_node_type::get_coor_s(R[k], nodes[k]);
 
  519            DIRECTOR_ESTIMATION::compute_normal(R, D);
 
  520            for (
int k = 0; k != nb_nodes; ++k) {
 
  521                if (!nodes_type_6_dof[k]
 
  522                    && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
  523                    if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
  526                    std::copy(R[k] + 3, R[k] + 6, D[k]);
 
  529            b2linalg::Vector<double> N(nb_nodes);
 
  530            SHAPE::eval_h(&internal_coor[0], N);
 
  531            property->laminate(N, e_begin, e_end);
 
  533            std::fill_n(D[0], 3 * nb_nodes, 0.0);
 
  541                z = face_id == 5 ? e_begin : e_end;
 
  542                map_node = SHAPE::face_node[face_id - 5];
 
  545                map_node = SHAPE::face_node[face_id - 7];
 
  547            nb_r_nodes = nb_nodes;
 
  549            map_node = SHAPE::edge_node[face_id - 1];
 
  550            nb_r_nodes = edge_shape::num_nodes;
 
  551            z = 0.5 * ((e_end - e_begin) * internal_coor[1] + e_end + e_begin);
 
  554        if (!geom.is_null()) {
 
  555            b2linalg::Vector<double> shape(nb_r_nodes);
 
  557                SHAPE::eval_h(&internal_coor[0], shape);
 
  559                edge_shape::eval_h(&internal_coor[0], shape);
 
  562            for (
size_t i = 0; i != 3; ++i) {
 
  564                for (
int k = 0; k != nb_r_nodes; ++k) {
 
  565                    size_t kk = map_node[k];
 
  566                    tmp += shape[k] * (R[kk][i] + z * D[kk][i]);
 
  571        if (!d_geom_d_icoor.is_null()) {
 
  573                b2linalg::Matrix<double> d_shape(2, nb_r_nodes);
 
  574                SHAPE::eval_h_derived(&internal_coor[0], d_shape);
 
  575                d_geom_d_icoor.resize(3, 2);
 
  576                for (
size_t i = 0; i != 3; ++i) {
 
  577                    for (
size_t j = 0; j != 2; ++j) {
 
  579                        for (
int k = 0; k != nb_r_nodes; ++k) {
 
  580                            size_t kk = map_node[k];
 
  581                            tmp += d_shape(j, k) * (R[kk][i] + z * D[kk][i]);
 
  583                        d_geom_d_icoor(i, j) = tmp;
 
  587                b2linalg::Matrix<double> d_shape(1, nb_r_nodes);
 
  588                edge_shape::eval_h_derived(&internal_coor[0], d_shape);
 
  589                b2linalg::Vector<double> shape(nb_r_nodes);
 
  590                edge_shape::eval_h(&internal_coor[0], shape);
 
  591                d_geom_d_icoor.resize(3, 2);
 
  592                const TP z_scale = 0.5 * (e_end - e_begin);
 
  593                for (
size_t i = 0; i != 3; ++i) {
 
  596                    for (
int k = 0; k != nb_r_nodes; ++k) {
 
  597                        size_t kk = map_node[k];
 
  598                        tmp0 += d_shape(0, k) * R[kk][i];
 
  599                        tmp1 += shape[k] * D[kk][i];
 
  601                    d_geom_d_icoor(i, 0) = tmp0;
 
  602                    d_geom_d_icoor(i, 1) = tmp1 * z_scale;
 
  608    void face_field_value(
 
  609          const int face_id, 
const std::string& field_name,
 
  610          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
 
  611          const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof, 
const double time,
 
  612          b2linalg::Vector<TP, b2linalg::Vdense>& value,
 
  613          b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_value_d_icoor,
 
  614          b2linalg::Index& dof_numbering, b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_value_d_dof,
 
  615          b2linalg::Index& d_value_d_dof_dep_col) {
 
  618        typedef typename SHAPE::edge_shape_type_0 edge_shape;
 
  619        TP R[nb_nodes][6] = {};
 
  620        TP D[nb_nodes][3] = {};
 
  623        for (
int k = 0; k != nb_nodes; ++k) {
 
  624            if (nodes_type_6_dof[k]) {
 
  625                dof6_node_type::get_coor_s(R[k], nodes[k]);
 
  627                dof5_node_type::get_coor_s(R[k], nodes[k]);
 
  630        DIRECTOR_ESTIMATION::compute_normal(R, D);
 
  631        for (
int k = 0; k != nb_nodes; ++k) {
 
  632            if (!nodes_type_6_dof[k]
 
  633                && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
  634                if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
  637                std::copy(R[k] + 3, R[k] + 6, D[k]);
 
  646                b2linalg::Vector<double> N(nb_nodes);
 
  647                SHAPE::eval_h(&internal_coor[0], N);
 
  648                property->laminate(N, e_begin, e_end);
 
  650                    z = 0.5 * ((e_end - e_begin) * internal_coor[1] + e_end + e_begin);
 
  651                    z_scale = 0.5 * (e_end - e_begin);
 
  653                    z = face_id == 5 ? e_begin : e_end;
 
  663            nb_nodes_f = edge_shape::num_nodes;
 
  664            map_node = SHAPE::edge_node[face_id - 1];
 
  667            nb_nodes_f = SHAPE::num_nodes;
 
  669                map_node = SHAPE::face_node[face_id - 5];
 
  671                map_node = SHAPE::face_node[face_id - 7];
 
  678                for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
  679                    std::copy(R[map_node[i]], R[map_node[i] + 1], tmp[i]);
 
  681                std::copy(tmp[0], tmp[edge_shape::num_nodes], R[0]);
 
  685                for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
  686                    std::copy(D[map_node[i]], D[map_node[i] + 1], tmp[i]);
 
  688                std::copy(tmp[0], tmp[edge_shape::num_nodes], D[0]);
 
  692        TP Dr[nb_nodes][2][3];
 
  695        for (
int k = 0; k != nb_nodes_f; ++k) {
 
  696            if (!nodes_type_6_dof[map_node[k]]) {
 
  697                const TP e2[3] = {0, 1, 0};
 
  700                    const TP e2[3] = {1, 0, 0};
 
  710        TP u_dof[nb_nodes][3];
 
  713        TP w_dof[nb_nodes][3];
 
  718        std::vector<size_t> dnvec(nb_nodes_f * 6 + 4);
 
  719        size_t* dn{dnvec.data()};
 
  723        std::fill_n(u_dof[0], nb_nodes_f * 3, 0);
 
  724        if (dof.size2() == 0) {
 
  725            std::fill_n(w_dof[0], nb_nodes_f * 3, 0);
 
  726            if (incompatible_mode) { std::fill_n(inc_dof, 4, 0); }
 
  727            for (
int k = 0; k != nb_nodes_f; ++k) {
 
  728                const int kk = map_node[k];
 
  729                if (nodes_type_6_dof[kk]) {
 
  730                    ptr_dn = dof6_node_type::get_global_dof_numbering_s(ptr_dn, nodes[kk]);
 
  731                    if (face_id > 6) { ptr_dn -= 3; }
 
  733                    ptr_dn = dof5_node_type::get_global_dof_numbering_s(ptr_dn, nodes[kk]);
 
  734                    if (face_id > 6) { ptr_dn -= 2; }
 
  738            for (
int k = 0; k != nb_nodes_f; ++k) {
 
  739                const int kk = map_node[k];
 
  740                if (nodes_type_6_dof[kk]) {
 
  741                    dof6_node_type::get_global_dof_numbering_s(ptr_dn, nodes[kk]);
 
  743                        for (
int i = 0; i != 3; ++i, ++ptr_dn) { u_dof[k][i] = dof(*ptr_dn, 0); }
 
  744                        for (
int i = 0; i != 3; ++i, ++ptr_dn) { w_dof[k][i] = dof(*ptr_dn, 0); }
 
  748                    if (face_id > 6) { ptr_dn -= 3; }
 
  750                    dof5_node_type::get_global_dof_numbering_s(ptr_dn, nodes[kk]);
 
  752                        for (
int i = 0; i != 3; ++i, ++ptr_dn) { u_dof[k][i] = dof(*ptr_dn, 0); }
 
  753                        const TP d0 = dof(*ptr_dn, 0);
 
  755                        const TP d1 = dof(*ptr_dn, 0);
 
  757                        w_dof[k][0] = Dr[k][0][0] * d0 + Dr[k][1][0] * d1;
 
  758                        w_dof[k][1] = Dr[k][0][1] * d0 + Dr[k][1][1] * d1;
 
  759                        w_dof[k][2] = Dr[k][0][2] * d0 + Dr[k][1][2] * d1;
 
  763                    if (face_id > 6) { ptr_dn -= 2; }
 
  780        TP T[nb_nodes][3][3];
 
  784        for (
int k = 0; k != nb_nodes_f; ++k) {
 
  786            const TP* w = w_dof[k];
 
  787            const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
  788            norm_w = std::sqrt(ww);
 
  791                const TP w0_2 = w[0] * w[0];
 
  792                const TP w1_2 = w[1] * w[1];
 
  793                const TP w2_2 = w[2] * w[2];
 
  794                O2[0] = -w2_2 - w1_2;
 
  797                O2[3] = -w2_2 - w0_2;
 
  799                O2[5] = -w1_2 - w0_2;
 
  809                s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
  810                c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
  811                cc = 1.0 / 6 - ww / 120 + w4 / 5040 - w6 / 362880 + w8 / 39916800;
 
  813                s = std::sin(norm_w) / norm_w;
 
  814                c = std::sin(norm_w * 0.5) / norm_w;
 
  816                cc = (norm_w - std::sin(norm_w)) / (ww * norm_w);
 
  820                const TP* 
const Dk = D[k];
 
  821                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
  822                        + (s * w[1] + c * O2[2]) * Dk[2];
 
  823                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
  824                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
  825                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
  826                        + (1 + c * O2[5]) * Dk[2];
 
  829            if (d_value_d_dof.is_null()) { 
continue; }
 
  834            if (nodes_type_6_dof[map_node[k]]) {
 
  835                HT3k[0][0] = 1 + c * O2[0];
 
  836                HT3k[1][0] = s * -w[2] + c * O2[1];
 
  837                HT3k[2][0] = s * w[1] + c * O2[2];
 
  838                HT3k[0][1] = s * w[2] + c * O2[1];
 
  839                HT3k[1][1] = 1 + c * O2[3];
 
  840                HT3k[2][1] = s * -w[0] + c * O2[4];
 
  841                HT3k[0][2] = s * -w[1] + c * O2[2];
 
  842                HT3k[1][2] = s * w[0] + c * O2[4];
 
  843                HT3k[2][2] = 1 + c * O2[5];
 
  846                    const TP* 
const Dk = Dr[k][0];
 
  847                    HT3k[0][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
  848                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
  849                    HT3k[0][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
  850                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
  851                    HT3k[0][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
  852                                 + (1 + c * O2[5]) * Dk[2];
 
  855                    const TP* 
const Dk = Dr[k][1];
 
  856                    HT3k[1][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
  857                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
  858                    HT3k[1][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
  859                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
  860                    HT3k[1][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
  861                                 + (1 + c * O2[5]) * Dk[2];
 
  865                const TP* 
const dk = d[k];
 
  866                T[k][0][0] = dk[2] * HT3k[0][1] - dk[1] * HT3k[0][2];
 
  867                T[k][0][1] = -dk[2] * HT3k[0][0] + dk[0] * HT3k[0][2];
 
  868                T[k][0][2] = dk[1] * HT3k[0][0] - dk[0] * HT3k[0][1];
 
  869                T[k][1][0] = dk[2] * HT3k[1][1] - dk[1] * HT3k[1][2];
 
  870                T[k][1][1] = -dk[2] * HT3k[1][0] + dk[0] * HT3k[1][2];
 
  871                T[k][1][2] = dk[1] * HT3k[1][0] - dk[0] * HT3k[1][1];
 
  872                if (nodes_type_6_dof[map_node[k]]) {
 
  873                    T[k][2][0] = dk[2] * HT3k[2][1] - dk[1] * HT3k[2][2];
 
  874                    T[k][2][1] = -dk[2] * HT3k[2][0] + dk[0] * HT3k[2][2];
 
  875                    T[k][2][2] = dk[1] * HT3k[2][0] - dk[0] * HT3k[2][1];
 
  880        if (!value.is_null()) {
 
  881            b2linalg::Vector<double> shape(nb_nodes_f);
 
  883                edge_shape::eval_h(&internal_coor[0], shape);
 
  885                SHAPE::eval_h(&internal_coor[0], shape);
 
  888            for (
size_t i = 0; i != 3; ++i) {
 
  890                for (
int k = 0; k != nb_nodes_f; ++k) {
 
  891                    tmp += shape[k] * (u_dof[k][i] + (d[k][i] - D[k][i]) * z);
 
  896        if (!d_value_d_icoor.is_null()) {
 
  898                b2linalg::Matrix<double> d_shape(2, nb_nodes_f);
 
  899                SHAPE::eval_h_derived(&internal_coor[0], d_shape);
 
  900                d_value_d_icoor.resize(3, 2);
 
  901                for (
size_t i = 0; i != 3; ++i) {
 
  902                    for (
size_t j = 0; j != 2; ++j) {
 
  904                        for (
int k = 0; k != nb_nodes_f; ++k) {
 
  905                            tmp += d_shape(j, k) * (u_dof[k][i] + (d[k][i] - D[k][i]) * z);
 
  907                        d_value_d_icoor(i, j) = tmp;
 
  911                b2linalg::Matrix<double> d_shape(1, nb_nodes_f);
 
  912                edge_shape::eval_h_derived(&internal_coor[0], d_shape);
 
  913                b2linalg::Vector<double> shape(nb_nodes_f);
 
  914                edge_shape::eval_h(&internal_coor[0], shape);
 
  915                d_value_d_icoor.resize(3, 2);
 
  916                for (
size_t i = 0; i != 3; ++i) {
 
  919                    for (
int k = 0; k != nb_nodes_f; ++k) {
 
  920                        tmp0 += d_shape(0, k) * u_dof[k][i];
 
  921                        tmp1 += shape[k] * (d[k][i] - D[k][i]);
 
  923                    d_value_d_icoor(i, 0) = tmp0;
 
  924                    d_value_d_icoor(i, 1) = tmp1 * z_scale;
 
  928        if (!d_value_d_dof.is_null()) {
 
  929            dof_numbering.assign(dn, ptr_dn);
 
  930            b2linalg::Vector<double> shape(nb_nodes_f);
 
  932                edge_shape::eval_h(&internal_coor[0], shape);
 
  934                SHAPE::eval_h(&internal_coor[0], shape);
 
  936            d_value_d_dof.resize(3, dof_numbering.size());
 
  937            d_value_d_dof.set_zero();
 
  939            for (
int k = 0; k != nb_nodes_f; ++k) {
 
  940                for (
size_t i = 0; i != 3; ++i) { d_value_d_dof(i, i_ptr + i) = shape[k]; }
 
  943                    size_t rs = nodes_type_6_dof[map_node[k]] ? 3 : 2;
 
  944                    for (
size_t j = 0; j != rs; ++j, ++i_ptr) {
 
  945                        for (
size_t i = 0; i != 3; ++i) {
 
  946                            d_value_d_dof(i, i_ptr) = shape[k] * z * T[k][j][i];
 
  951            if (!d_value_d_dof_dep_col.is_null() && face_id < 7) {
 
  954                    for (
int k = 0; k != nb_nodes_f; ++k) {
 
  955                        if (nodes_type_6_dof[map_node[k]]) {
 
  958                            for (
int i = 1; i != 3; ++i) {
 
  959                                if (std::abs(R[k][i]) > std::abs(R[k][i_max])) { i_max = i; }
 
  961                            d_value_d_dof_dep_col.push_back(i_ptr + i_max);
 
  968                    for (
int k = 0; k != nb_nodes_f; ++k) {
 
  969                        const int i_max = nodes_type_6_dof[map_node[k]] ? 3 : 2;
 
  971                        for (
int i = 0; i != i_max; ++i, ++i_ptr) {
 
  972                            d_value_d_dof_dep_col.push_back(i_ptr);
 
  981          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
 
  982          b2linalg::Vector<TP>& geom, b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_geom_d_icoor) {
 
  986        TP D[nb_nodes][3] = {};
 
  989        TP R[nb_nodes][6] = {};
 
  991        for (
int k = 0; k != nb_nodes_ip; ++k) {
 
  992            if (nodes_type_6_dof[k]) {
 
  993                dof6_node_type::get_coor_s(R[k], nodes[k]);
 
  995                dof5_node_type::get_coor_s(R[k], nodes[k]);
 
  999        for (
int k = nb_nodes_ip; k != nb_nodes; ++k) {
 
 1000            dof2_node_type::get_coor_s(R[k], nodes[k]);
 
 1004        DIRECTOR_ESTIMATION::compute_normal(R, D);
 
 1005        for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 1006            if (!nodes_type_6_dof[k]
 
 1007                && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
 1008                if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
 1011                std::copy(R[k] + 3, R[k] + 6, D[k]);
 
 1017            b2linalg::Vector<double> N(nb_nodes);
 
 1018            SHAPE::eval_h(&internal_coor[0], N);
 
 1020            property->laminate(N, e_begin, e_end);
 
 1021            z = 0.5 * ((e_end - e_begin) * internal_coor[2] + e_end + e_begin);
 
 1022            dz = 0.5 * (e_end - e_begin);
 
 1025        if (!geom.is_null()) {
 
 1027            b2linalg::Vector<double> shape(nb_nodes);
 
 1028            SHAPE::eval_h(&internal_coor[0], shape);
 
 1029            for (
size_t i = 0; i != 3; ++i) {
 
 1031                for (
int k = 0; k != nb_nodes; ++k) { tmp += shape[k] * (R[k][i] + D[k][i] * z); }
 
 1035        if (!d_geom_d_icoor.is_null()) {
 
 1036            b2linalg::Vector<double> shape(nb_nodes);
 
 1037            SHAPE::eval_h(&internal_coor[0], shape);
 
 1038            b2linalg::Matrix<double> d_shape(2, nb_nodes);
 
 1039            SHAPE::eval_h_derived(&internal_coor[0], d_shape);
 
 1040            d_geom_d_icoor.resize(3, 3);
 
 1041            for (
size_t i = 0; i != 3; ++i) {
 
 1042                for (
size_t j = 0; j != 2; ++j) {
 
 1044                    for (
int k = 0; k != nb_nodes; ++k) {
 
 1045                        tmp += d_shape(j, k) * (R[k][i] + D[k][i] * z);
 
 1047                    d_geom_d_icoor(i, j) = tmp;
 
 1051                    for (
int k = 0; k != nb_nodes; ++k) { tmp += shape[k] * D[k][i] * dz; }
 
 1052                    d_geom_d_icoor(i, 2) = tmp;
 
 1058    bool body_field_linear_on_dof(
const std::string& field_name)
 override { 
return false; }
 
 1060    void body_field_value(
 
 1061          const std::string& field_name,
 
 1062          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
 
 1063          const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof, 
const double time,
 
 1064          b2linalg::Vector<TP, b2linalg::Vdense>& value,
 
 1065          b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_value_d_icoor,
 
 1066          b2linalg::Index& dof_numbering,
 
 1067          b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_value_d_dof) {
 
 1071        TP D[nb_nodes][3] = {};
 
 1074        TP Dr[nb_nodes][2][3] = {};
 
 1077        TP R[nb_nodes][6] = {};
 
 1079        int number_of_dof = incompatible_mode ? 4 : 0;
 
 1081        for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 1082            if (nodes_type_6_dof[k]) {
 
 1084                dof6_node_type::get_coor_s(R[k], nodes[k]);
 
 1087                dof5_node_type::get_coor_s(R[k], nodes[k]);
 
 1091        for (
int k = nb_nodes_ip; k != nb_nodes; ++k) {
 
 1093            dof2_node_type::get_coor_s(R[k], nodes[k]);
 
 1097        DIRECTOR_ESTIMATION::compute_normal(R, D);
 
 1098        for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 1099            if (!nodes_type_6_dof[k]
 
 1100                && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
 1101                if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
 1104                std::copy(R[k] + 3, R[k] + 6, D[k]);
 
 1109        for (
int k = 0; k != nb_nodes; ++k) {
 
 1110            if (!nodes_type_6_dof[k]) {
 
 1111                const TP e2[3] = {0, 1, 0};
 
 1114                    const TP e2[3] = {1, 0, 0};
 
 1124        TP u_dof[nb_nodes][3];
 
 1127        TP w_dof[nb_nodes][3];
 
 1130        if (dof.size2() == 0) {
 
 1131            std::fill_n(u_dof[0], nb_nodes * 3, 0);
 
 1132            std::fill_n(w_dof[0], nb_nodes * 3, 0);
 
 1134            for (
int k = 0; k != nb_nodes; ++k) {
 
 1136                nodes[k]->get_global_dof_numbering(dn);
 
 1137                for (
int i = 0; i != 3; ++i) { u_dof[k][i] = dof(dn[i], 0); }
 
 1138                if (nodes_type_6_dof[k]) {
 
 1139                    for (
int i = 0; i != 3; ++i) { w_dof[k][i] = dof(dn[3 + i], 0); }
 
 1141                    const TP d0 = dof(dn[3], 0);
 
 1142                    const TP d1 = dof(dn[4], 0);
 
 1143                    w_dof[k][0] = Dr[k][0][0] * d0 + Dr[k][1][0] * d1;
 
 1144                    w_dof[k][1] = Dr[k][0][1] * d0 + Dr[k][1][1] * d1;
 
 1145                    w_dof[k][2] = Dr[k][0][2] * d0 + Dr[k][1][2] * d1;
 
 1155        TP HT3[nb_nodes][3][3];
 
 1156        TP T[nb_nodes][3][3];
 
 1157        TP norm_w[nb_nodes];
 
 1161        for (
int k = 0; k != nb_nodes; ++k) {
 
 1162            const TP* w = w_dof[k];
 
 1163            const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
 1164            norm_w[k] = std::sqrt(ww);
 
 1167                const TP w0_2 = w[0] * w[0];
 
 1168                const TP w1_2 = w[1] * w[1];
 
 1169                const TP w2_2 = w[2] * w[2];
 
 1170                O2[0] = -w2_2 - w1_2;
 
 1171                O2[1] = w[0] * w[1];
 
 1172                O2[2] = w[0] * w[2];
 
 1173                O2[3] = -w2_2 - w0_2;
 
 1174                O2[4] = w[1] * w[2];
 
 1175                O2[5] = -w1_2 - w0_2;
 
 1185                s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
 1186                c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
 1187                cc = 1.0 / 6 - ww / 120 + w4 / 5040 - w6 / 362880 + w8 / 39916800;
 
 1189                s = std::sin(norm_w[k]) / norm_w[k];
 
 1190                c = std::sin(norm_w[k] * 0.5) / norm_w[k];
 
 1192                cc = (norm_w[k] - std::sin(norm_w[k])) / (ww * norm_w[k]);
 
 1195                TP* 
const dk = d[k];
 
 1196                const TP* 
const Dk = D[k];
 
 1197                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 1198                        + (s * w[1] + c * O2[2]) * Dk[2];
 
 1199                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 1200                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
 1201                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 1202                        + (1 + c * O2[5]) * Dk[2];
 
 1205            if (!d_value_d_dof.is_null()) {
 
 1209                typedef TP HT3k_t[3][3];
 
 1210                HT3k_t& HT3k = HT3[k];
 
 1211                if (nodes_type_6_dof[k]) {
 
 1212                    HT3k[0][0] = 1 + c * O2[0];
 
 1213                    HT3k[1][0] = s * -w[2] + c * O2[1];
 
 1214                    HT3k[2][0] = s * w[1] + c * O2[2];
 
 1215                    HT3k[0][1] = s * w[2] + c * O2[1];
 
 1216                    HT3k[1][1] = 1 + c * O2[3];
 
 1217                    HT3k[2][1] = s * -w[0] + c * O2[4];
 
 1218                    HT3k[0][2] = s * -w[1] + c * O2[2];
 
 1219                    HT3k[1][2] = s * w[0] + c * O2[4];
 
 1220                    HT3k[2][2] = 1 + c * O2[5];
 
 1223                        const TP* 
const Dk = Dr[k][0];
 
 1224                        HT3k[0][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 1225                                     + (s * w[1] + c * O2[2]) * Dk[2];
 
 1226                        HT3k[0][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 1227                                     + (s * -w[0] + c * O2[4]) * Dk[2];
 
 1228                        HT3k[0][2] = (s * -w[1] + c * O2[2]) * Dk[0]
 
 1229                                     + (s * w[0] + c * O2[4]) * Dk[1] + (1 + c * O2[5]) * Dk[2];
 
 1232                        const TP* 
const Dk = Dr[k][1];
 
 1233                        HT3k[1][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 1234                                     + (s * w[1] + c * O2[2]) * Dk[2];
 
 1235                        HT3k[1][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 1236                                     + (s * -w[0] + c * O2[4]) * Dk[2];
 
 1237                        HT3k[1][2] = (s * -w[1] + c * O2[2]) * Dk[0]
 
 1238                                     + (s * w[0] + c * O2[4]) * Dk[1] + (1 + c * O2[5]) * Dk[2];
 
 1242                    const TP* 
const dk = d[k];
 
 1243                    T[k][0][0] = dk[2] * HT3k[0][1] - dk[1] * HT3k[0][2];
 
 1244                    T[k][0][1] = -dk[2] * HT3k[0][0] + dk[0] * HT3k[0][2];
 
 1245                    T[k][0][2] = dk[1] * HT3k[0][0] - dk[0] * HT3k[0][1];
 
 1246                    T[k][1][0] = dk[2] * HT3k[1][1] - dk[1] * HT3k[1][2];
 
 1247                    T[k][1][1] = -dk[2] * HT3k[1][0] + dk[0] * HT3k[1][2];
 
 1248                    T[k][1][2] = dk[1] * HT3k[1][0] - dk[0] * HT3k[1][1];
 
 1249                    if (nodes_type_6_dof[k]) {
 
 1250                        T[k][2][0] = dk[2] * HT3k[2][1] - dk[1] * HT3k[2][2];
 
 1251                        T[k][2][1] = -dk[2] * HT3k[2][0] + dk[0] * HT3k[2][2];
 
 1252                        T[k][2][2] = dk[1] * HT3k[2][0] - dk[0] * HT3k[2][1];
 
 1260            b2linalg::Vector<double> N(nb_nodes);
 
 1261            SHAPE::eval_h(&internal_coor[0], N);
 
 1263            property->laminate(N, e_begin, e_end);
 
 1264            z = 0.5 * ((e_end - e_begin) * internal_coor[2] + e_end + e_begin);
 
 1265            dz = 0.5 * (e_end - e_begin);
 
 1268        if (!value.is_null()) {
 
 1269            b2linalg::Vector<double> shape(nb_nodes);
 
 1270            SHAPE::eval_h(&internal_coor[0], shape);
 
 1272            for (
size_t i = 0; i != 3; ++i) {
 
 1274                for (
int k = 0; k != nb_nodes; ++k) {
 
 1275                    tmp += shape[k] * (u_dof[k][i] + (d[k][i] - D[k][i]) * z);
 
 1280        if (!d_value_d_icoor.is_null()) {
 
 1281            b2linalg::Vector<double> shape(nb_nodes);
 
 1282            SHAPE::eval_h(&internal_coor[0], shape);
 
 1283            b2linalg::Matrix<double> d_shape(2, nb_nodes);
 
 1284            SHAPE::eval_h_derived(&internal_coor[0], d_shape);
 
 1285            d_value_d_icoor.resize(3, 3);
 
 1286            for (
size_t i = 0; i != 3; ++i) {
 
 1287                for (
size_t j = 0; j != 2; ++j) {
 
 1289                    for (
int k = 0; k != nb_nodes; ++k) {
 
 1290                        tmp += d_shape(j, k) * (u_dof[k][i] + (d[k][i] - D[k][i]) * z);
 
 1292                    d_value_d_icoor(i, j) = tmp;
 
 1296                    for (
int k = 0; k != nb_nodes; ++k) {
 
 1297                        tmp += shape[k] * (d[k][i] - D[k][i]) * dz;
 
 1299                    d_value_d_icoor(i, 2) = tmp;
 
 1303        if (!d_value_d_dof.is_null()) {
 
 1304            get_dof_numbering(dof_numbering);
 
 1305            b2linalg::Vector<double> shape(nb_nodes);
 
 1306            SHAPE::eval_h(&internal_coor[0], shape);
 
 1307            d_value_d_dof.resize(3, dof_numbering.size());
 
 1308            d_value_d_dof.set_zero();
 
 1310            for (
int k = 0; k != nb_nodes; ++k) {
 
 1311                for (
size_t i = 0; i != 3; ++i, ++i_ptr) { d_value_d_dof(i, i_ptr) = shape[k]; }
 
 1312                size_t rs = nodes_type_6_dof[k] ? 3 : 2;
 
 1313                for (
size_t j = 0; j != rs; ++j, ++i_ptr) {
 
 1314                    for (
size_t i = 0; i != 3; ++i) {
 
 1315                        d_value_d_dof(i, i_ptr) = shape[k] * z * T[k][j][i];
 
 1322    void field_edge_integration(
 
 1323          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dof,
 
 1324          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdot,
 
 1325          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdotdot, 
const Field<TP>& f,
 
 1326          int edge, b2linalg::Index& dof_numbering,
 
 1327          b2linalg::Vector<TP, b2linalg::Vdense>& discretised_field,
 
 1328          b2linalg::Matrix<TP, b2linalg::Mpacked>& d_discretised_field_d_dof);
 
 1330    void field_face_integration(
 
 1331          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dof,
 
 1332          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdot,
 
 1333          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdotdot, 
const Field<TP>& field,
 
 1334          int face, b2linalg::Index& dof_numbering,
 
 1335          b2linalg::Vector<TP, b2linalg::Vdense>& discretised_field,
 
 1336          b2linalg::Matrix<TP, b2linalg::Mpacked>& d_discretised_field_d_dof);
 
 1338    void field_volume_integration(
 
 1339          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dof,
 
 1340          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdot,
 
 1341          const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdotdot, 
const Field<TP>& field,
 
 1342          b2linalg::Index& dof_numbering, b2linalg::Vector<TP, b2linalg::Vdense>& discretised_field,
 
 1343          b2linalg::Matrix<TP, b2linalg::Mpacked>& d_discretised_field_d_dof);
 
 1348    static const int nb_nodes = SHAPE::num_nodes;
 
 1349    static const int nb_nodes_ip = SHAPE_IP::num_nodes;
 
 1352#ifdef TT_INTEGRATION_AT_GAUSS_POINT 
 1353#if TT_INTEGRATION_AT_GAUSS_POINT == 2 
 1354    static const int nb_thickness_gauss = 2;
 
 1355#elif TT_INTEGRATION_AT_GAUSS_POINT == 3 
 1356    static const int nb_thickness_gauss = 3;
 
 1359    static const int nb_thickness_gauss = 2;
 
 1362    static const int max_nb_dof = nb_nodes * 6 + (incompatible_mode ? 4 : 0);
 
 1365    size_t incompatible_dof_index;
 
 1368    Node* nodes[nb_nodes];
 
 1371    SolidMechanics25D<TP>* property;
 
 1372    TP non_structural_mass;
 
 1374    bool line_load_as_moment;
 
 1376    int get_number_of_all_dof()
 const {
 
 1377        return 5 * nb_nodes_ip + nodes_type_6_dof.count() + 2 * (nb_nodes - nb_nodes_ip)
 
 1378               + (incompatible_mode ? 4 : 0);
 
 1383    std::bitset<nb_nodes> nodes_type_6_dof;
 
 1386    typedef SolidMechanics25D<TP>* properties_list_t[nb_gauss][nb_thickness_gauss];
 
 1390        properties_list_t* properties_list;
 
 1398    struct OnePointShape {
 
 1406        double dN[2][nb_nodes];
 
 1409        double N_ip[nb_nodes_ip];
 
 1412        double dN_ip[2][nb_nodes_ip];
 
 1415    struct OnePoint : 
public OnePointShape, 
public MITC_TYING_POINT::TyingPointInterp {
 
 1422    static OnePoint gauss_point[nb_gauss];
 
 1425    static OnePointShape tying_point[MITC_TYING_POINT::nb_tying_point];
 
 1428    static OnePointShape central_point;
 
 1431    static void init_shape_point(
double r, 
double s, OnePointShape& ops) {
 
 1432        b2linalg::Vector<double, b2linalg::Vdense> v(nb_nodes);
 
 1433        b2linalg::Matrix<double, b2linalg::Mrectangle> m(2, nb_nodes);
 
 1436        SHAPE::eval_h(ops.IntCoor, v);
 
 1437        std::copy(v.begin(), v.end(), ops.N);
 
 1438        SHAPE::eval_h_derived(ops.IntCoor, m);
 
 1439        for (
size_t i = 0; i != m.size1(); ++i) {
 
 1440            for (
size_t j = 0; j != m.size2(); ++j) { ops.dN[i][j] = m(i, j); }
 
 1442        m.resize(2, nb_nodes_ip);
 
 1443        SHAPE_IP::eval_h(ops.IntCoor, v);
 
 1444        std::copy(v.begin(), v.end(), ops.N_ip);
 
 1445        SHAPE::eval_h_derived(ops.IntCoor, m);
 
 1446        for (
size_t i = 0; i != m.size1(); ++i) {
 
 1447            for (
size_t j = 0; j != m.size2(); ++j) { ops.dN_ip[i][j] = m(i, j); }
 
 1451    static void init_gauss_point(
double r, 
double s, OnePoint& op) {
 
 1452        init_shape_point(r, s, op);
 
 1453        MITC_TYING_POINT::get_interpolation(r, s, op);
 
 1456    static void init_all_tying_point() {
 
 1457        for (
int i = 0; i != MITC_TYING_POINT::nb_tying_point; ++i) {
 
 1459                  MITC_TYING_POINT::tying_point[i].r, MITC_TYING_POINT::tying_point[i].s,
 
 1464    static void init_all_gauss_point() {
 
 1465        INTEGRATION integration;
 
 1466        integration.set_num(nb_gauss);
 
 1467        if (nb_gauss != integration.get_num()) {
 
 1468            Exception() << 
"MITC shell element  did not find an integration " 
 1469                           "schema that contains integration points." 
 1473        for (
int l = 0; l != nb_gauss; ++l) {
 
 1474            integration.get_point(l, ip_id, gauss_point[l].weight);
 
 1475            init_gauss_point(ip_id[0], ip_id[1], gauss_point[l]);
 
 1479#ifdef CHECK_TYING_INTERPOLATION 
 1481    static void check_tying_interpolation() {
 
 1482        const double tol = 1e-10;
 
 1483        std::cerr << 
"Check element " 
 1484                  << (NAME_SUFFIX ? std::string(SHAPE::name) + 
".S.MITC" + NAME_SUFFIX
 
 1485                                  : std::string(SHAPE::name) + 
".S.MITC")
 
 1487        for (
int i = 0; i != MITC_TYING_POINT::nb_tying_point; ++i) {
 
 1488            typename MITC_TYING_POINT::TyingPointInterp interp;
 
 1489            MITC_TYING_POINT::get_interpolation(
 
 1490                  MITC_TYING_POINT::tying_point[i].r, MITC_TYING_POINT::tying_point[i].s, interp);
 
 1491            for (
int j = 0; j != MITC_TYING_POINT::nb_err_comp; ++j) {
 
 1492                if (interp.err[j].tp_pos == i) {
 
 1493                    if (std::abs(interp.err[j].weight - 1) > tol) {
 
 1494                        std::cerr << 
"err[" << j << 
"] at " << j << 
" = " << interp.err[j].weight
 
 1497                    for (
int k = 0; k != MITC_TYING_POINT::nb_err_comp; ++k) {
 
 1498                        if (k != j && std::abs(interp.err[k].weight) > tol) {
 
 1499                            std::cerr << 
"err[" << k << 
"] at " << j << 
" = " 
 1500                                      << interp.err[k].weight << std::endl;
 
 1506            for (
int j = 0; j != MITC_TYING_POINT::nb_ess_comp; ++j) {
 
 1507                if (interp.ess[j].tp_pos == i) {
 
 1508                    if (std::abs(interp.ess[j].weight - 1) > tol) {
 
 1509                        std::cerr << 
"ess[" << j << 
"] at " << j << 
" = " << interp.ess[j].weight
 
 1512                    for (
int k = 0; k != MITC_TYING_POINT::nb_ess_comp; ++k) {
 
 1513                        if (k != j && std::abs(interp.ess[k].weight) > tol) {
 
 1514                            std::cerr << 
"ess[" << k << 
"] at " << j << 
" = " 
 1515                                      << interp.ess[k].weight << std::endl;
 
 1521            for (
int j = 0; j != MITC_TYING_POINT::nb_ers_comp; ++j) {
 
 1522                if (interp.ers[j].tp_pos == i) {
 
 1523                    if (std::abs(interp.ers[j].weight - 1) > tol) {
 
 1524                        std::cerr << 
"ers[" << j << 
"] at " << j << 
" = " << interp.ers[j].weight
 
 1527                    for (
int k = 0; k != MITC_TYING_POINT::nb_ers_comp; ++k) {
 
 1528                        if (k != j && std::abs(interp.ers[k].weight) > tol) {
 
 1529                            std::cerr << 
"ers[" << k << 
"] at " << j << 
" = " 
 1530                                      << interp.ers[k].weight << std::endl;
 
 1536            for (
int j = 0; j != MITC_TYING_POINT::nb_ert_comp; ++j) {
 
 1537                if (interp.ert[j].tp_pos == i) {
 
 1538                    if (std::abs(interp.ert[j].weight - 1) > tol) {
 
 1539                        std::cerr << 
"ert[" << j << 
"] at " << j << 
" = " << interp.ert[j].weight
 
 1542                    for (
int k = 0; k != MITC_TYING_POINT::nb_ert_comp; ++k) {
 
 1543                        if (k != j && std::abs(interp.ert[k].weight) > tol) {
 
 1544                            std::cerr << 
"ert[" << k << 
"] at " << j << 
" = " 
 1545                                      << interp.ert[k].weight << std::endl;
 
 1551            for (
int j = 0; j != MITC_TYING_POINT::nb_est_comp; ++j) {
 
 1552                if (interp.est[j].tp_pos == i) {
 
 1553                    if (std::abs(interp.est[j].weight - 1) > tol) {
 
 1554                        std::cerr << 
"est[" << j << 
"] at " << j << 
" = " << interp.est[j].weight
 
 1557                    for (
int k = 0; k != MITC_TYING_POINT::nb_est_comp; ++k) {
 
 1558                        if (k != j && std::abs(interp.est[k].weight) > tol) {
 
 1559                            std::cerr << 
"est[" << k << 
"] at " << j << 
" = " 
 1560                                      << interp.est[k].weight << std::endl;
 
 1570    struct InitElemType {
 
 1575#ifdef CHECK_TYING_INTERPOLATION 
 1576            ShellMITC::check_tying_interpolation();
 
 1579            ShellMITC::init_all_tying_point();
 
 1580            ShellMITC::init_all_gauss_point();
 
 1582            init_shape_point(0, 0, central_point);
 
 1587    static const InitElemType init_elem_type;
 
 1591      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1592      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1593      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1595      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1596      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::nb_nodes;
 
 1599      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1600      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1601      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1603      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1604      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::OnePoint
 
 1606            TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1607            incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::gauss_point[nb_gauss];
 
 1610      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1611      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1612      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1614      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1615      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::OnePointShape
 
 1617            TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1618            incompatible_mode, mitc_gs, nb_gauss, INTEGRATION,
 
 1619            NAME_SUFFIX>::tying_point[MITC_TYING_POINT::nb_tying_point];
 
 1623      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1624      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1625      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1627      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1628      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::OnePointShape
 
 1630            TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1631            incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::central_point;
 
 1635      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1636      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1637      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1638const typename ShellMITC<
 
 1639      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1640      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::InitElemType
 
 1642            TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1643            incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::init_elem_type;
 
 1646      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1647      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1648      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1650      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1651      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::type_t
 
 1653            TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1654            incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::
 
 1655            type((NAME_SUFFIX ? std::string(SHAPE::name) + 
".S.MITC" + NAME_SUFFIX
 
 1656                              : std::string(SHAPE::name) + 
".S.MITC")
 
 1657                       + (incompatible_mode ? 
".E4" : 
"") + (mitc_gs ? 
".GS" : 
"")
 
 1658                       + (std::is_same<TP, 
b2000::csda<double>>::value ? 
".CSDA" : 
""),
 
 1659                 "", StringList(), element::module, &TypedElement<double>::type);
 
 1664      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 1665      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 1666      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 1667void b2000::ShellMITC<
 
 1668      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 1669      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::
 
 1671            Model& model, 
const bool linear, 
const EquilibriumSolution equilibrium_solution,
 
 1672            const double time, 
const double delta_time,
 
 1673            const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof,
 
 1674            GradientContainer* gradient_container, SolverHints* solver_hints,
 
 1675            b2linalg::Index& dof_numbering, b2linalg::Vector<TP, b2linalg::Vdense>& value,
 
 1676            const std::vector<bool>& d_value_d_dof_flags,
 
 1677            std::vector<b2linalg::Matrix<TP, b2linalg::Mpacked>>& d_value_d_dof,
 
 1678            b2linalg::Vector<TP, b2linalg::Vdense>& d_value_d_time) {
 
 1679#ifdef TT_INTEGRATION_AT_GAUSS_POINT 
 1680#if TT_INTEGRATION_AT_GAUSS_POINT == 2 
 1681    const double sqrt_inv_3 = std::sqrt(1 / 3.);
 
 1682    static const int nb_tt_integration = 2;
 
 1683    static const double n_omega[2] = {-sqrt_inv_3, sqrt_inv_3};
 
 1684    static const double n_weight[2] = {0.5, 0.5};
 
 1685#elif TT_INTEGRATION_AT_GAUSS_POINT == 3 
 1686    static const int nb_tt_integration = 3;
 
 1687    static const double n_omega[3] = {-1, 0, 1};
 
 1688    static const double n_weight[3] = {1. / 6, 4. / 6, 1. / 6};
 
 1691    static const int nb_tt_integration = 2;
 
 1692    static const double n_omega[2] = {-1, 1};
 
 1693    static const double n_weight[2] = {0.5, 0.5};
 
 1696    get_dof_numbering(dof_numbering);
 
 1697    if (!value.is_null()) {
 
 1698        value.resize(dof_numbering.size());
 
 1701    assert(d_value_d_dof.size() == d_value_d_dof_flags.size());
 
 1702    for (
size_t i = 0; i != d_value_d_dof.size(); ++i) {
 
 1703        if (d_value_d_dof_flags[i]) {
 
 1704            d_value_d_dof[i].resize(dof_numbering.size(), 
true);
 
 1705            d_value_d_dof[i].set_zero();
 
 1707            d_value_d_dof[i].resize(
size_t(0), 
true);
 
 1710    if (!d_value_d_time.is_null()) {
 
 1711        d_value_d_time.resize(dof_numbering.size());
 
 1712        d_value_d_time.set_zero();
 
 1716    TP D[nb_nodes][3] = {};
 
 1717    TP Dd[nb_nodes][3] = {};
 
 1720    TP Dr[nb_nodes][2][3] = {};
 
 1723    TP R[nb_nodes][6] = {};
 
 1725    const int number_of_dof = (int)dof_numbering.size();
 
 1729    for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 1730        if (nodes_type_6_dof[k]) {
 
 1731            dof6_node_type::get_coor_s(R[k], nodes[k]);
 
 1733            dof5_node_type::get_coor_s(R[k], nodes[k]);
 
 1736    for (
int k = nb_nodes_ip; k != nb_nodes; ++k) { dof2_node_type::get_coor_s(R[k], nodes[k]); }
 
 1738    bool drill_stabilized_node[nb_nodes] = {};
 
 1741    DIRECTOR_ESTIMATION::compute_normal(R, D);
 
 1742    for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 1743        if (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5) {
 
 1744            if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
 1747            std::copy(R[k] + 3, R[k] + 6, D[k]);
 
 1748            drill_stabilized_node[k] = nodes_type_6_dof[k];
 
 1752            drill_stabilized_node[k] = 
false;
 
 1757    for (
int k = 0; k != nb_nodes; ++k) {
 
 1758        if (!nodes_type_6_dof[k]) {
 
 1759            const TP e2[3] = {0, 1, 0};
 
 1762                const TP e2[3] = {1, 0, 0};
 
 1772    TP u_dof[nb_nodes][3];
 
 1775    TP w_dof[nb_nodes][3];
 
 1781    if (dof.size2() == 0 || dof.size1() == 0) {
 
 1782        std::fill_n(u_dof[0], nb_nodes * 3, 0);
 
 1783        std::fill_n(w_dof[0], nb_nodes * 3, 0);
 
 1784        if (incompatible_mode) { std::fill_n(inc_dof, 4, 0); }
 
 1787        for (
int k = 0; k != nb_nodes; ++k) {
 
 1788            const TP* alpha = &dof[0][dof_numbering[o]];
 
 1789            if (nodes_type_6_dof[k]) {
 
 1790                std::copy(alpha, alpha + 3, u_dof[k]);
 
 1791                std::copy(alpha + 3, alpha + 6, w_dof[k]);
 
 1793            } 
else if (SHAPE_IP_OOP_SAME || k < nb_nodes_ip) {
 
 1794                std::copy(alpha, alpha + 3, u_dof[k]);
 
 1795                w_dof[k][0] = Dr[k][0][0] * alpha[3] + Dr[k][1][0] * alpha[4];
 
 1796                w_dof[k][1] = Dr[k][0][1] * alpha[3] + Dr[k][1][1] * alpha[4];
 
 1797                w_dof[k][2] = Dr[k][0][2] * alpha[3] + Dr[k][1][2] * alpha[4];
 
 1800                w_dof[k][0] = Dr[k][0][0] * alpha[0] + Dr[k][1][0] * alpha[1];
 
 1801                w_dof[k][1] = Dr[k][0][1] * alpha[0] + Dr[k][1][1] * alpha[1];
 
 1802                w_dof[k][2] = Dr[k][0][2] * alpha[0] + Dr[k][1][2] * alpha[1];
 
 1806        if (incompatible_mode) {
 
 1807            const TP* alpha = &dof[0][dof_numbering[o]];
 
 1808            std::copy(alpha, alpha + 4, inc_dof);
 
 1817    TP HT3[nb_nodes][3][3];
 
 1818    TP T[nb_nodes][3][3];
 
 1819    TP norm_w[nb_nodes];
 
 1823    for (
int k = 0; k != nb_nodes; ++k) {
 
 1824        const TP* w = w_dof[k];
 
 1825        const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
 1826        norm_w[k] = std::sqrt(ww);
 
 1829            const TP w0_2 = w[0] * w[0];
 
 1830            const TP w1_2 = w[1] * w[1];
 
 1831            const TP w2_2 = w[2] * w[2];
 
 1832            O2[0] = -w2_2 - w1_2;
 
 1833            O2[1] = w[0] * w[1];
 
 1834            O2[2] = w[0] * w[2];
 
 1835            O2[3] = -w2_2 - w0_2;
 
 1836            O2[4] = w[1] * w[2];
 
 1837            O2[5] = -w1_2 - w0_2;
 
 1847            s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
 1848            c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
 1849            cc = 1.0 / 6 - ww / 120 + w4 / 5040 - w6 / 362880 + w8 / 39916800;
 
 1851            s = std::sin(norm_w[k]) / norm_w[k];
 
 1852            c = std::sin(norm_w[k] * 0.5) / norm_w[k];
 
 1854            cc = (norm_w[k] - std::sin(norm_w[k])) / (ww * norm_w[k]);
 
 1857            TP* 
const dk = d[k];
 
 1858            const TP* 
const Dk = D[k];
 
 1860                TP* 
const Ddk = Dd[k];
 
 1861                Ddk[0] = Dk[2] * w_dof[k][1] - Dk[1] * w_dof[k][2];
 
 1862                Ddk[1] = -Dk[2] * w_dof[k][0] + Dk[0] * w_dof[k][2];
 
 1863                Ddk[2] = Dk[1] * w_dof[k][0] - Dk[0] * w_dof[k][1];
 
 1864                dk[0] = Ddk[0] + Dk[0];
 
 1865                dk[1] = Ddk[1] + Dk[1];
 
 1866                dk[2] = Ddk[2] + Dk[2];
 
 1868                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 1869                        + (s * w[1] + c * O2[2]) * Dk[2];
 
 1870                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 1871                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
 1872                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 1873                        + (1 + c * O2[5]) * Dk[2];
 
 1874                TP* 
const Ddk = Dd[k];
 
 1875                Ddk[0] = dk[0] - Dk[0];
 
 1876                Ddk[1] = dk[1] - Dk[1];
 
 1877                Ddk[2] = dk[2] - Dk[2];
 
 1882            const TP* 
const dk = D[k];
 
 1883            if (nodes_type_6_dof[k]) {
 
 1885                T[k][0][1] = -dk[2];
 
 1889                T[k][1][2] = -dk[0];
 
 1890                T[k][2][0] = -dk[1];
 
 1894                T[k][0][0] = dk[2] * Dr[k][0][1] - dk[1] * Dr[k][0][2];
 
 1895                T[k][0][1] = -dk[2] * Dr[k][0][0] + dk[0] * Dr[k][0][2];
 
 1896                T[k][0][2] = dk[1] * Dr[k][0][0] - dk[0] * Dr[k][0][1];
 
 1897                T[k][1][0] = dk[2] * Dr[k][1][1] - dk[1] * Dr[k][1][2];
 
 1898                T[k][1][1] = -dk[2] * Dr[k][1][0] + dk[0] * Dr[k][1][2];
 
 1899                T[k][1][2] = dk[1] * Dr[k][1][0] - dk[0] * Dr[k][1][1];
 
 1905            typedef TP HT3k_t[3][3];
 
 1906            HT3k_t& HT3k = HT3[k];
 
 1907            if (nodes_type_6_dof[k]) {
 
 1908                HT3k[0][0] = 1 + c * O2[0];
 
 1909                HT3k[1][0] = s * -w[2] + c * O2[1];
 
 1910                HT3k[2][0] = s * w[1] + c * O2[2];
 
 1911                HT3k[0][1] = s * w[2] + c * O2[1];
 
 1912                HT3k[1][1] = 1 + c * O2[3];
 
 1913                HT3k[2][1] = s * -w[0] + c * O2[4];
 
 1914                HT3k[0][2] = s * -w[1] + c * O2[2];
 
 1915                HT3k[1][2] = s * w[0] + c * O2[4];
 
 1916                HT3k[2][2] = 1 + c * O2[5];
 
 1919                    const TP* 
const Dk = Dr[k][0];
 
 1920                    HT3k[0][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 1921                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
 1922                    HT3k[0][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 1923                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
 1924                    HT3k[0][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 1925                                 + (1 + c * O2[5]) * Dk[2];
 
 1928                    const TP* 
const Dk = Dr[k][1];
 
 1929                    HT3k[1][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 1930                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
 1931                    HT3k[1][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 1932                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
 1933                    HT3k[1][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 1934                                 + (1 + c * O2[5]) * Dk[2];
 
 1938                const TP* 
const dk = d[k];
 
 1939                T[k][0][0] = dk[2] * HT3k[0][1] - dk[1] * HT3k[0][2];
 
 1940                T[k][0][1] = -dk[2] * HT3k[0][0] + dk[0] * HT3k[0][2];
 
 1941                T[k][0][2] = dk[1] * HT3k[0][0] - dk[0] * HT3k[0][1];
 
 1942                T[k][1][0] = dk[2] * HT3k[1][1] - dk[1] * HT3k[1][2];
 
 1943                T[k][1][1] = -dk[2] * HT3k[1][0] + dk[0] * HT3k[1][2];
 
 1944                T[k][1][2] = dk[1] * HT3k[1][0] - dk[0] * HT3k[1][1];
 
 1945                if (nodes_type_6_dof[k]) {
 
 1946                    T[k][2][0] = dk[2] * HT3k[2][1] - dk[1] * HT3k[2][2];
 
 1947                    T[k][2][1] = -dk[2] * HT3k[2][0] + dk[0] * HT3k[2][2];
 
 1948                    T[k][2][2] = dk[1] * HT3k[2][0] - dk[0] * HT3k[2][1];
 
 1963        TP contravariant[mitc_gs ? 3 : 0][3];
 
 1965    TyingPoint tp_value[MITC_TYING_POINT::nb_tying_point];
 
 1968    for (
int l = 0; l != MITC_TYING_POINT::nb_tying_point; ++l) {
 
 1969        TyingPoint& tp = tp_value[l];
 
 1970        const OnePointShape& N = tying_point[l];
 
 1981        for (
int k = 0; k != nb_nodes; ++k) {
 
 1982            for (
int j = 0; j != 3; ++j) {
 
 1983                Ar[j] += N.dN[0][k] * R[k][j];
 
 1984                As[j] += N.dN[1][k] * R[k][j];
 
 1985                At[j] += N.N[k] * D[k][j];
 
 1986                Atr[j] += N.dN[0][k] * D[k][j];
 
 1987                Ats[j] += N.dN[1][k] * D[k][j];
 
 1988                if (SHAPE_IP_OOP_SAME) {
 
 1989                    ur[j] += N.dN[0][k] * u_dof[k][j];
 
 1990                    us[j] += N.dN[1][k] * u_dof[k][j];
 
 1992                ut[j] += N.N[k] * Dd[k][j];
 
 1993                utr[j] += N.dN[0][k] * Dd[k][j];
 
 1994                uts[j] += N.dN[1][k] * Dd[k][j];
 
 1997        if (!SHAPE_IP_OOP_SAME) {
 
 1998            for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 1999                for (
int j = 0; j != 3; ++j) {
 
 2000                    ur[j] += N.dN_ip[0][k] * u_dof[k][j];
 
 2001                    us[j] += N.dN_ip[1][k] * u_dof[k][j];
 
 2006            for (
int j = 0; j != 3; ++j) {
 
 2007                tp.ars[0][j] = Ar[j];
 
 2008                tp.ars[1][j] = As[j];
 
 2010                tp.atrs[0][j] = Atr[j];
 
 2011                tp.atrs[1][j] = Ats[j];
 
 2014            for (
int j = 0; j != 3; ++j) {
 
 2015                tp.ars[0][j] = Ar[j] + ur[j];
 
 2016                tp.ars[1][j] = As[j] + us[j];
 
 2017                tp.at[j] = At[j] + ut[j];
 
 2018                tp.atrs[0][j] = Atr[j] + utr[j];
 
 2019                tp.atrs[1][j] = Ats[j] + uts[j];
 
 2024            for (
int j = 0; j != 3; ++j) {
 
 2031        if (MITC_TYING_POINT::tying_point[l].comp_ip || mitc_gs) {
 
 2033            tp.Ers0[0] = Ar[0] * ur[0] + Ar[1] * ur[1] + Ar[2] * ur[2];
 
 2035            tp.Ers0[1] = As[0] * us[0] + As[1] * us[1] + As[2] * us[2];
 
 2037            tp.Ers0[2] = Ar[0] * us[0] + Ar[1] * us[1] + Ar[2] * us[2] + As[0] * ur[0]
 
 2038                         + As[1] * ur[1] + As[2] * ur[2];
 
 2040            tp.Ers1[0] = Ar[0] * utr[0] + Ar[1] * utr[1] + Ar[2] * utr[2] + ur[0] * Atr[0]
 
 2041                         + ur[1] * Atr[1] + ur[2] * Atr[2];
 
 2043            tp.Ers1[1] = As[0] * uts[0] + As[1] * uts[1] + As[2] * uts[2] + us[0] * Ats[0]
 
 2044                         + us[1] * Ats[1] + us[2] * Ats[2];
 
 2046            tp.Ers1[2] = Ar[0] * uts[0] + Ar[1] * uts[1] + Ar[2] * uts[2] + ur[0] * Ats[0]
 
 2047                         + ur[1] * Ats[1] + ur[2] * Ats[2] + As[0] * utr[0] + As[1] * utr[1]
 
 2048                         + As[2] * utr[2] + us[0] * Atr[0] + us[1] * Atr[1] + us[2] * Atr[2];
 
 2050                tp.Ers0[0] += 0.5 * (ur[0] * ur[0] + ur[1] * ur[1] + ur[2] * ur[2]);
 
 2051                tp.Ers0[1] += 0.5 * (us[0] * us[0] + us[1] * us[1] + us[2] * us[2]);
 
 2052                tp.Ers0[2] += ur[0] * us[0] + ur[1] * us[1] + ur[2] * us[2];
 
 2053                tp.Ers1[0] += ur[0] * utr[0] + ur[1] * utr[1] + ur[2] * utr[2];
 
 2054                tp.Ers1[1] += us[0] * uts[0] + us[1] * uts[1] + us[2] * uts[2];
 
 2055                tp.Ers1[2] += ur[0] * uts[0] + ur[1] * uts[1] + ur[2] * uts[2] + us[0] * utr[0]
 
 2056                              + us[1] * utr[1] + us[2] * utr[2];
 
 2059            if (incompatible_mode) {
 
 2061                const TP a00 = -2 * N.IntCoor[0] * inc_dof[0];
 
 2062                const TP a10 = -2 * N.IntCoor[0] * inc_dof[2];
 
 2063                const TP a01 = -2 * N.IntCoor[1] * inc_dof[1];
 
 2064                const TP a11 = -2 * N.IntCoor[1] * inc_dof[3];
 
 2068                    tp.Ers0[2] += a01 + a10;
 
 2070                    tp.Ers0[0] += a00 + 0.5 * (a00 * a00 + a10 * a10);
 
 2071                    tp.Ers0[1] += a11 + 0.5 * (a11 * a11 + a01 * a01);
 
 2072                    tp.Ers0[2] += a01 + a10 + a00 * a01 + a10 * a11;
 
 2076        if (MITC_TYING_POINT::tying_point[l].comp_oop || mitc_gs) {
 
 2078            tp.Erst[0] = At[0] * ur[0] + At[1] * ur[1] + At[2] * ur[2] + ut[0] * Ar[0]
 
 2079                         + ut[1] * Ar[1] + ut[2] * Ar[2];
 
 2081            tp.Erst[1] = At[0] * us[0] + At[1] * us[1] + At[2] * us[2] + ut[0] * As[0]
 
 2082                         + ut[1] * As[1] + ut[2] * As[2];
 
 2084                tp.Erst[0] += ut[0] * ur[0] + ut[1] * ur[1] + ut[2] * ur[2];
 
 2085                tp.Erst[1] += ut[0] * us[0] + ut[1] * us[1] + ut[2] * us[2];
 
 2097    TP Bl[max_nb_dof][nb_gauss][8];
 
 2099    const bool linear_mat = 
property->linear() == SolidMechanics::yes;
 
 2100    const bool path_dependent_mat = 
property->path_dependent() && !linear;
 
 2101    bool need_to_compute_C_linear = 
false;
 
 2102    const bool compute_value_from_dofdotdot = !value.is_null() && !dof.is_null()
 
 2103                                              && d_value_d_dof_flags.size() > 2
 
 2104                                              && d_value_d_dof_flags[2];
 
 2105    const bool d_value_d_dofdotdot_null =
 
 2106          d_value_d_dof_flags.size() <= 2 || !d_value_d_dof_flags[2];
 
 2107    const bool compute_stress =
 
 2108          !value.is_null() || (!linear && d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0])
 
 2109          || (equilibrium_solution && path_dependent_mat);
 
 2110    const bool property_shell_stress = 
property->has_shell_stress_interface();
 
 2115    std::fill_n(h[0], nb_nodes * 3, 0);
 
 2117    TP C_tmp[nb_gauss][36];
 
 2119    if (linear_mat && !linear && !property_shell_stress) {
 
 2121            C_linear = 
new TP[nb_gauss][36];
 
 2122            std::fill_n(C_linear[0], 36 * nb_gauss, 0);
 
 2123            need_to_compute_C_linear = 
true;
 
 2127        if (d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]) {
 
 2129            need_to_compute_C_linear = 
true;
 
 2130            std::fill_n(C[0], 36 * nb_gauss, 0);
 
 2132        if (!value.is_null() || compute_stress) { std::fill_n(S[0], 8 * nb_gauss, 0); }
 
 2150    TP tp_trans[mitc_gs ? nb_gauss : 0][MITC_TYING_POINT::nb_tying_point][5][5];
 
 2153    for (
int l = 0; l != nb_gauss; ++l) {
 
 2154        const OnePoint& N = gauss_point[l];
 
 2175        for (
int k = 0; k != nb_nodes; ++k) {
 
 2176            for (
int j = 0; j != 3; ++j) {
 
 2177                Ar[j] += N.dN[0][k] * R[k][j];
 
 2178                As[j] += N.dN[1][k] * R[k][j];
 
 2179                At[j] += N.N[k] * D[k][j];
 
 2180                Atr[j] += N.dN[0][k] * D[k][j];
 
 2181                Ats[j] += N.dN[1][k] * D[k][j];
 
 2182                if (SHAPE_IP_OOP_SAME) {
 
 2183                    ur[j] += N.dN[0][k] * u_dof[k][j];
 
 2184                    us[j] += N.dN[1][k] * u_dof[k][j];
 
 2186                ut[j] += N.N[k] * Dd[k][j];
 
 2187                utr[j] += N.dN[0][k] * Dd[k][j];
 
 2188                uts[j] += N.dN[1][k] * Dd[k][j];
 
 2191        if (!SHAPE_IP_OOP_SAME) {
 
 2192            for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 2193                for (
int j = 0; j != 3; ++j) {
 
 2194                    ur[j] += N.dN_ip[0][k] * u_dof[k][j];
 
 2195                    us[j] += N.dN_ip[1][k] * u_dof[k][j];
 
 2202            for (
int j = 0; j != 3; ++j) {
 
 2207            for (
int tk = 0; tk != MITC_TYING_POINT::nb_tying_point; ++tk) {
 
 2210                TP* C = tp_trans[l][tk][0];
 
 2211                C[0] = trans[0][0] * trans[0][0];
 
 2212                C[1] = trans[0][1] * trans[0][1];
 
 2213                C[2] = trans[0][0] * trans[0][1];
 
 2214                C[3] = trans[0][2] * trans[0][0];
 
 2215                C[4] = trans[0][1] * trans[0][2];
 
 2217                C = tp_trans[l][tk][1];
 
 2218                C[0] = trans[1][0] * trans[1][0];
 
 2219                C[1] = trans[1][1] * trans[1][1];
 
 2220                C[2] = trans[1][0] * trans[1][1];
 
 2221                C[3] = trans[1][2] * trans[1][0];
 
 2222                C[4] = trans[1][1] * trans[1][2];
 
 2224                C = tp_trans[l][tk][2];
 
 2225                C[0] = 2 * trans[0][0] * trans[1][0];
 
 2226                C[1] = 2 * trans[0][1] * trans[1][1];
 
 2227                C[2] = (trans[0][0] * trans[1][1] + trans[0][1] * trans[1][0]);
 
 2228                C[3] = (trans[0][2] * trans[1][0] + trans[0][0] * trans[1][2]);
 
 2229                C[4] = (trans[0][1] * trans[1][2] + trans[0][2] * trans[1][1]);
 
 2231                C = tp_trans[l][tk][3];
 
 2232                C[0] = 2 * trans[2][0] * trans[0][0];
 
 2233                C[1] = 2 * trans[2][1] * trans[0][1];
 
 2234                C[2] = (trans[2][0] * trans[0][1] + trans[2][1] * trans[0][0]);
 
 2235                C[3] = (trans[2][2] * trans[0][0] + trans[2][0] * trans[0][2]);
 
 2236                C[4] = (trans[2][1] * trans[0][2] + trans[2][2] * trans[0][1]);
 
 2238                C = tp_trans[l][tk][4];
 
 2239                C[0] = 2 * trans[1][0] * trans[2][0];
 
 2240                C[1] = 2 * trans[1][1] * trans[2][1];
 
 2241                C[2] = (trans[1][0] * trans[2][1] + trans[1][1] * trans[2][0]);
 
 2242                C[3] = (trans[1][2] * trans[2][0] + trans[1][0] * trans[2][2]);
 
 2243                C[4] = (trans[1][1] * trans[2][2] + trans[1][2] * trans[2][1]);
 
 2248        if (gradient_container || compute_stress) {
 
 2249            TP* 
const El = E[l];
 
 2251            if (MITC_TYING_POINT::nb_err_comp == 0) {
 
 2252                El[0] = Ar[0] * ur[0] + Ar[1] * ur[1] + Ar[2] * ur[2];
 
 2253                El[3] = Ar[0] * utr[0] + Ar[1] * utr[1] + Ar[2] * utr[2] + ur[0] * Atr[0]
 
 2254                        + ur[1] * Atr[1] + ur[2] * Atr[2];
 
 2258                for (
int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
 
 2259                    const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
 
 2261                        El[0] += interp.weight * tp_value[interp.tp_pos].Ers0[interp.e_pos];
 
 2262                        El[3] += interp.weight * tp_value[interp.tp_pos].Ers1[interp.e_pos];
 
 2264                        const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 2265                        El[0] += interp.weight
 
 2266                                 * (tpt[0] * tp_value[interp.tp_pos].Ers0[0]
 
 2267                                    + tpt[1] * tp_value[interp.tp_pos].Ers0[1]
 
 2268                                    + tpt[2] * tp_value[interp.tp_pos].Ers0[2]
 
 2269                                    + tpt[3] * tp_value[interp.tp_pos].Erst[0]
 
 2270                                    + tpt[4] * tp_value[interp.tp_pos].Erst[1]);
 
 2271                        El[3] += interp.weight
 
 2272                                 * (tpt[0] * tp_value[interp.tp_pos].Ers1[0]
 
 2273                                    + tpt[1] * tp_value[interp.tp_pos].Ers1[1]
 
 2274                                    + tpt[2] * tp_value[interp.tp_pos].Ers1[2]);
 
 2279            if (MITC_TYING_POINT::nb_ess_comp == 0) {
 
 2280                El[1] = As[0] * us[0] + As[1] * us[1] + As[2] * us[2];
 
 2281                El[4] = As[0] * uts[0] + As[1] * uts[1] + As[2] * uts[2] + us[0] * Ats[0]
 
 2282                        + us[1] * Ats[1] + us[2] * Ats[2];
 
 2286                for (
int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
 
 2287                    const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
 
 2289                        El[1] += interp.weight * tp_value[interp.tp_pos].Ers0[interp.e_pos];
 
 2290                        El[4] += interp.weight * tp_value[interp.tp_pos].Ers1[interp.e_pos];
 
 2292                        const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 2293                        El[1] += interp.weight
 
 2294                                 * (tpt[0] * tp_value[interp.tp_pos].Ers0[0]
 
 2295                                    + tpt[1] * tp_value[interp.tp_pos].Ers0[1]
 
 2296                                    + tpt[2] * tp_value[interp.tp_pos].Ers0[2]
 
 2297                                    + tpt[3] * tp_value[interp.tp_pos].Erst[0]
 
 2298                                    + tpt[4] * tp_value[interp.tp_pos].Erst[1]);
 
 2299                        El[4] += interp.weight
 
 2300                                 * (tpt[0] * tp_value[interp.tp_pos].Ers1[0]
 
 2301                                    + tpt[1] * tp_value[interp.tp_pos].Ers1[1]
 
 2302                                    + tpt[2] * tp_value[interp.tp_pos].Ers1[2]);
 
 2307            if (MITC_TYING_POINT::nb_ers_comp == 0) {
 
 2308                El[2] = Ar[0] * us[0] + Ar[1] * us[1] + Ar[2] * us[2] + As[0] * ur[0]
 
 2309                        + As[1] * ur[1] + As[2] * ur[2];
 
 2310                El[5] = Ar[0] * uts[0] + Ar[1] * uts[1] + Ar[2] * uts[2] + ur[0] * Ats[0]
 
 2311                        + ur[1] * Ats[1] + ur[2] * Ats[2] + As[0] * utr[0] + As[1] * utr[1]
 
 2312                        + As[2] * utr[2] + us[0] * Atr[0] + us[1] * Atr[1] + us[2] * Atr[2];
 
 2316                for (
int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
 
 2317                    const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
 
 2319                        El[2] += interp.weight * tp_value[interp.tp_pos].Ers0[interp.e_pos];
 
 2320                        El[5] += interp.weight * tp_value[interp.tp_pos].Ers1[interp.e_pos];
 
 2322                        const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 2323                        El[2] += interp.weight
 
 2324                                 * (tpt[0] * tp_value[interp.tp_pos].Ers0[0]
 
 2325                                    + tpt[1] * tp_value[interp.tp_pos].Ers0[1]
 
 2326                                    + tpt[2] * tp_value[interp.tp_pos].Ers0[2]
 
 2327                                    + tpt[3] * tp_value[interp.tp_pos].Erst[0]
 
 2328                                    + tpt[4] * tp_value[interp.tp_pos].Erst[1]);
 
 2329                        El[5] += interp.weight
 
 2330                                 * (tpt[0] * tp_value[interp.tp_pos].Ers1[0]
 
 2331                                    + tpt[1] * tp_value[interp.tp_pos].Ers1[1]
 
 2332                                    + tpt[2] * tp_value[interp.tp_pos].Ers1[2]);
 
 2337            if (MITC_TYING_POINT::nb_ert_comp == 0) {
 
 2338                El[6] = At[0] * ur[0] + At[1] * ur[1] + At[2] * ur[2] + ut[0] * Ar[0]
 
 2339                        + ut[1] * Ar[1] + ut[2] * Ar[2];
 
 2342                for (
int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
 
 2343                    const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
 
 2345                        El[6] += interp.weight * tp_value[interp.tp_pos].Erst[interp.e_pos];
 
 2347                        const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 2348                        El[6] += interp.weight
 
 2349                                 * (tpt[0] * tp_value[interp.tp_pos].Ers0[0]
 
 2350                                    + tpt[1] * tp_value[interp.tp_pos].Ers0[1]
 
 2351                                    + tpt[2] * tp_value[interp.tp_pos].Ers0[2]
 
 2352                                    + tpt[3] * tp_value[interp.tp_pos].Erst[0]
 
 2353                                    + tpt[4] * tp_value[interp.tp_pos].Erst[1]);
 
 2358            if (MITC_TYING_POINT::nb_est_comp == 0) {
 
 2359                El[7] = At[0] * us[0] + At[1] * us[1] + At[2] * us[2] + ut[0] * As[0]
 
 2360                        + ut[1] * As[1] + ut[2] * As[2];
 
 2363                for (
int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
 
 2364                    const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
 
 2366                        El[7] += interp.weight * tp_value[interp.tp_pos].Erst[interp.e_pos];
 
 2368                        const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 2369                        El[7] += interp.weight
 
 2370                                 * (tpt[0] * tp_value[interp.tp_pos].Ers0[0]
 
 2371                                    + tpt[1] * tp_value[interp.tp_pos].Ers0[1]
 
 2372                                    + tpt[2] * tp_value[interp.tp_pos].Ers0[2]
 
 2373                                    + tpt[3] * tp_value[interp.tp_pos].Erst[0]
 
 2374                                    + tpt[4] * tp_value[interp.tp_pos].Erst[1]);
 
 2379                if (MITC_TYING_POINT::nb_err_comp == 0) {
 
 2380                    El[0] += 0.5 * (ur[0] * ur[0] + ur[1] * ur[1] + ur[2] * ur[2]);
 
 2381                    El[3] += ur[0] * utr[0] + ur[1] * utr[1] + ur[2] * utr[2];
 
 2383                if (MITC_TYING_POINT::nb_ess_comp == 0) {
 
 2384                    El[1] += 0.5 * (us[0] * us[0] + us[1] * us[1] + us[2] * us[2]);
 
 2385                    El[4] += us[0] * uts[0] + us[1] * uts[1] + us[2] * uts[2];
 
 2387                if (MITC_TYING_POINT::nb_ers_comp == 0) {
 
 2388                    El[2] += ur[0] * us[0] + ur[1] * us[1] + ur[2] * us[2];
 
 2389                    El[5] += ur[0] * uts[0] + ur[1] * uts[1] + ur[2] * uts[2] + us[0] * utr[0]
 
 2390                             + us[1] * utr[1] + us[2] * utr[2];
 
 2392                if (MITC_TYING_POINT::nb_ert_comp == 0) {
 
 2393                    El[6] += ut[0] * ur[0] + ut[1] * ur[1] + ut[2] * ur[2];
 
 2395                if (MITC_TYING_POINT::nb_est_comp == 0) {
 
 2396                    El[7] += ut[0] * us[0] + ut[1] * us[1] + ut[2] * us[2];
 
 2399            if (incompatible_mode) {
 
 2401                const TP a00 = -2 * N.IntCoor[0] * inc_dof[0];
 
 2402                const TP a10 = -2 * N.IntCoor[0] * inc_dof[2];
 
 2403                const TP a01 = -2 * N.IntCoor[1] * inc_dof[1];
 
 2404                const TP a11 = -2 * N.IntCoor[1] * inc_dof[3];
 
 2410                    El[0] += a00 + 0.5 * (a00 * a00 + a10 * a10);
 
 2411                    El[1] += a11 + 0.5 * (a11 * a11 + a01 * a01);
 
 2412                    El[2] += a01 + a10 + a00 * a01 + a10 * a11;
 
 2419        for (
int i = 0; i != nb_nodes; ++i) {
 
 2420            coor[0] += N.N[i] * R[i][0];
 
 2421            coor[1] += N.N[i] * R[i][1];
 
 2422            coor[2] += N.N[i] * R[i][2];
 
 2425        if (gradient_container || (compute_stress && !(linear_mat && C != 0))
 
 2426            || need_to_compute_C_linear) {
 
 2427            if (property_shell_stress) {
 
 2430                for (
int j = 0; j != 3; ++j) {
 
 2446                if (gradient_container) {
 
 2448                    const GradientContainer::InternalElementPosition ip = {
 
 2449                          N.IntCoor[0], N.IntCoor[1], 0, 0};
 
 2450                    gradient_container->set_current_position(ip, w);
 
 2451                    static const GradientContainer::FieldDescription coor_descr = {
 
 2452                          "COOR_IP", 
"x.y.z",   
"Coordinate of the point", 3, 3, 1, 3,
 
 2454                    if (gradient_container->compute_field_value(coor_descr.name)) {
 
 2455                        gradient_container->set_field_value(coor_descr, coor);
 
 2460                    static const GradientContainer::FieldDescription disp_descr = {
 
 2461                          "DISP_IP", 
"dx.dy.dz", 
"Displacement at the point", 3, 3, 1, 3,
 
 2463                    if (gradient_container->compute_field_value(disp_descr.name)) {
 
 2465                        for (
int i = 0; i != nb_nodes; ++i) {
 
 2466                            disp[0] += N.N[i] * u_dof[i][0];
 
 2467                            disp[1] += N.N[i] * u_dof[i][1];
 
 2468                            disp[2] += N.N[i] * u_dof[i][2];
 
 2470                        gradient_container->set_field_value(disp_descr, disp);
 
 2472                    static const GradientContainer::FieldDescription vdescr = {
 
 2473                          "VOLUME", 
"V",           
"Integration point volume", 3, 1, -1, -1,
 
 2474                          false,    
typeid(double)};
 
 2475                    const double w_double = w;
 
 2476                    if (gradient_container->compute_field_value(vdescr.name)) {
 
 2477                        gradient_container->set_field_value(vdescr, &w_double);
 
 2481                const double el_coordinates[3] = {N.IntCoor[0], N.IntCoor[1], 0};
 
 2482                assert(!need_to_compute_C_linear || C != 0);
 
 2483                (path_dependent_mat ? properties_list[0][l][0] : property)
 
 2484                      ->get_static_nonlinear_shell_stress(
 
 2485                            &model, time, delta_time, equilibrium_solution, gradient_container,
 
 2486                            solver_hints, 
this, el_coordinates,
 
 2487                            b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N),
 
 2493                            (compute_stress ? S[l] : 0),            
 
 2494                            (need_to_compute_C_linear ? C[l] : 0),  
 
 2498                if (compute_stress) {
 
 2499                    for (
int i = 0; i != 8; ++i) { S[l][i] *= w; }
 
 2501                if (need_to_compute_C_linear) {
 
 2502                    for (
int i = 0; i != 36; ++i) { C[l][i] *= w; }
 
 2508                const int number_of_layer = 
property->number_of_layer();
 
 2510                      b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin,
 
 2513                for (
int layer = 0; layer != number_of_layer; ++layer) {
 
 2517                          b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
 
 2519                    for (
int lt = 0; lt != nb_tt_integration; ++lt) {
 
 2520                        TP omega = l_begin + 0.5 * (1 + n_omega[lt]) * (l_end - l_begin);
 
 2521                        const double eps = 1e-10;
 
 2522                        if (std::abs(omega + 1) < eps) {
 
 2524                        } 
else if (std::abs(omega) < eps) {
 
 2526                        } 
else if (std::abs(omega - 1) < eps) {
 
 2529                        const bool write_gradients = number_of_layer == 1 || n_omega[lt] == 0
 
 2531                                                     || (layer + 1 == number_of_layer);
 
 2535                        for (
int j = 0; j != 3; ++j) {
 
 2539                            u[0][j] = ur[j] + omega * utr[j];
 
 2540                            u[1][j] = us[j] + omega * uts[j];
 
 2542                            G[0][j] = Ar[j] + omega * Atr[j];
 
 2543                            G[1][j] = As[j] + omega * Ats[j];
 
 2548                        TP* 
const El = E[l];
 
 2549                        const TP EE[6] = {El[0] + omega * El[3],
 
 2550                                          El[1] + omega * El[4],
 
 2552                                          El[2] + omega * El[5],
 
 2555                        const double thickness_interpolation[2] = {
 
 2556                              (omega - e_begin) / (e_end - e_begin),
 
 2557                              (e_end - omega) / (e_end - e_begin)};
 
 2558                        const double el_coordinates[3] = {
 
 2559                              N.IntCoor[0], N.IntCoor[1], double(omega / (e_end - e_begin) * 2)};
 
 2560                        TP w = n_weight[lt] * (l_end - l_begin) * det_G * gauss_point[l].weight;
 
 2565                        if (gradient_container && write_gradients) {
 
 2567                            const GradientContainer::InternalElementPosition ip = {
 
 2568                                  N.IntCoor[0], N.IntCoor[1], omega / (e_end - e_begin) * 2, layer};
 
 2569                            gradient_container->set_current_position(ip, w);
 
 2570                            static const GradientContainer::FieldDescription coor_descr = {
 
 2571                                  "COOR_IP", 
"x.y.z",   
"Coordinate of the point", 3, 3, 1, 3,
 
 2573                            if (gradient_container->compute_field_value(coor_descr.name)) {
 
 2575                                for (
int i = 0; i != nb_nodes; ++i) {
 
 2576                                    coor[0] += N.N[i] * (R[i][0] + omega * D[i][0]);
 
 2577                                    coor[1] += N.N[i] * (R[i][1] + omega * D[i][1]);
 
 2578                                    coor[2] += N.N[i] * (R[i][2] + omega * D[i][2]);
 
 2580                                gradient_container->set_field_value(coor_descr, coor);
 
 2585                            static const GradientContainer::FieldDescription disp_descr = {
 
 2586                                  "DISP_IP", 
"dx.dy.dz", 
"Displacement at the point", 3, 3, 1, 3,
 
 2588                            if (gradient_container->compute_field_value(disp_descr.name)) {
 
 2590                                for (
int i = 0; i != nb_nodes; ++i) {
 
 2591                                    disp[0] += N.N[i] * (u_dof[i][0] + omega * (d[i][0] - D[i][0]));
 
 2592                                    disp[1] += N.N[i] * (u_dof[i][1] + omega * (d[i][1] - D[i][1]));
 
 2593                                    disp[2] += N.N[i] * (u_dof[i][2] + omega * (d[i][2] - D[i][2]));
 
 2595                                gradient_container->set_field_value(disp_descr, disp);
 
 2597                            static const GradientContainer::FieldDescription vdescr = {
 
 2598                                  "VOLUME", 
"V",           
"Integration point volume", 3, 1, -1, -1,
 
 2599                                  false,    
typeid(double)};
 
 2600                            const double w_double = w;
 
 2601                            if (gradient_container->compute_field_value(vdescr.name)) {
 
 2602                                gradient_container->set_field_value(vdescr, &w_double);
 
 2608                        if (need_to_compute_C_linear) {
 
 2615                        if (compute_stress && !(linear_mat && C != 0)) {
 
 2621                        (path_dependent_mat ? properties_list[layer][l][lt] : property)
 
 2623                                    &model, linear, equilibrium_solution, time, delta_time,
 
 2624                                    (write_gradients ? gradient_container : nullptr), solver_hints,
 
 2625                                    this, el_coordinates, layer,
 
 2626                                    b2linalg::Vector<double, b2linalg::Vdense_constref>(
 
 2628                                    thickness_interpolation,
 
 2651                                const TP ow = omega * w;
 
 2652                                Sl[3] += ow * Ss[0];
 
 2653                                Sl[4] += ow * Ss[1];
 
 2654                                Sl[5] += ow * Ss[3];
 
 2666                                Cl[15] += w * Cc[15];
 
 2668                                const TP ow = omega * w;
 
 2669                                Cl[3] += ow * Cc[0];
 
 2670                                Cl[4] += ow * Cc[1];
 
 2671                                Cl[5] += ow * Cc[3];
 
 2672                                Cl[10] += ow * Cc[1];
 
 2673                                Cl[11] += ow * Cc[6];
 
 2674                                Cl[12] += ow * Cc[8];
 
 2675                                Cl[16] += ow * Cc[3];
 
 2676                                Cl[17] += ow * Cc[8];
 
 2677                                Cl[18] += ow * Cc[15];
 
 2681                                Cl[13] += w * Cc[10];
 
 2682                                Cl[14] += w * Cc[9];
 
 2683                                Cl[19] += w * Cc[17];
 
 2684                                Cl[20] += w * Cc[16];
 
 2686                                const TP oow = ow * omega;
 
 2687                                Cl[21] += oow * Cc[0];
 
 2688                                Cl[22] += oow * Cc[1];
 
 2689                                Cl[23] += oow * Cc[3];
 
 2690                                Cl[26] += oow * Cc[6];
 
 2691                                Cl[27] += oow * Cc[8];
 
 2692                                Cl[30] += oow * Cc[15];
 
 2694                                Cl[24] += ow * Cc[5];
 
 2695                                Cl[25] += ow * Cc[4];
 
 2696                                Cl[28] += ow * Cc[10];
 
 2697                                Cl[29] += ow * Cc[9];
 
 2698                                Cl[31] += ow * Cc[17];
 
 2699                                Cl[32] += ow * Cc[16];
 
 2701                                Cl[33] += w * Cc[20];
 
 2702                                Cl[34] += w * Cc[19];
 
 2703                                Cl[35] += w * Cc[18];
 
 2712        TP J_tim[incompatible_mode ? 2 : 0][incompatible_mode ? 2 : 0];
 
 2713        if (incompatible_mode) {}
 
 2716        if ((d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]) || !value.is_null()) {
 
 2717            if (compute_stress && linear_mat && C != 0) {
 
 2719                blas::spmv(
'L', 8, 1.0, C[l], E[l], 1, 0.0, S[l], 1);
 
 2723                for (
int j = 0; j != 3; ++j) {
 
 2731                for (
int j = 0; j != 3; ++j) {
 
 2732                    ar[j] = Ar[j] + ur[j];
 
 2733                    as[j] = As[j] + us[j];
 
 2734                    at[j] = At[j] + ut[j];
 
 2735                    atr[j] = Atr[j] + utr[j];
 
 2736                    ats[j] = Ats[j] + uts[j];
 
 2742            for (
int k = 0; k != nb_nodes; ++k) {
 
 2743                if (SHAPE_IP_OOP_SAME || k < nb_nodes_ip) {
 
 2744                    for (
int j = 0; j != 3; ++j, ++dofi) {
 
 2745                        TP* 
const Blp = Bl[dofi][l];
 
 2746                        if (MITC_TYING_POINT::nb_err_comp == 0) {
 
 2747                            Blp[0] = N.dN_ip[0][k] * ar[j];
 
 2748                            Blp[3] = N.dN_ip[0][k] * atr[j];
 
 2752                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
 
 2753                                const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
 
 2755                                    if (interp.e_pos != 2) {
 
 2758                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 2759                                              * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 2762                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 2763                                              * tp_value[interp.tp_pos].atrs[interp.e_pos][j];
 
 2765                                        Blp[0] += interp.weight
 
 2766                                                  * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 2767                                                           * tp_value[interp.tp_pos].ars[1][j]
 
 2768                                                     + tying_point[interp.tp_pos].dN_ip[1][k]
 
 2769                                                             * tp_value[interp.tp_pos].ars[0][j]);
 
 2770                                        Blp[3] += interp.weight
 
 2771                                                  * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 2772                                                           * tp_value[interp.tp_pos].atrs[1][j]
 
 2773                                                     + tying_point[interp.tp_pos].dN_ip[1][k]
 
 2774                                                             * tp_value[interp.tp_pos].atrs[0][j]);
 
 2777                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 2780                                          * (tp_value[interp.tp_pos].ars[0][j]
 
 2782                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2784                                                              * tying_point[interp.tp_pos]
 
 2786                                             + tp_value[interp.tp_pos].ars[1][j]
 
 2788                                                              * tying_point[interp.tp_pos]
 
 2791                                                                * tying_point[interp.tp_pos]
 
 2793                                             + tp_value[interp.tp_pos].at[j]
 
 2794                                                     * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
 
 2796                                                                * tying_point[interp.tp_pos]
 
 2800                                          * (tp_value[interp.tp_pos].atrs[0][j]
 
 2802                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2804                                                              * tying_point[interp.tp_pos]
 
 2806                                             + tp_value[interp.tp_pos].atrs[1][j]
 
 2808                                                              * tying_point[interp.tp_pos]
 
 2811                                                                * tying_point[interp.tp_pos]
 
 2816                        if (MITC_TYING_POINT::nb_ess_comp == 0) {
 
 2817                            Blp[1] = N.dN_ip[1][k] * as[j];
 
 2818                            Blp[4] = N.dN_ip[1][k] * ats[j];
 
 2822                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
 
 2823                                const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
 
 2825                                    if (interp.e_pos != 2) {
 
 2828                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 2829                                              * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 2832                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 2833                                              * tp_value[interp.tp_pos].atrs[interp.e_pos][j];
 
 2835                                        Blp[1] += interp.weight
 
 2836                                                  * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 2837                                                           * tp_value[interp.tp_pos].ars[1][j]
 
 2838                                                     + tying_point[interp.tp_pos].dN_ip[1][k]
 
 2839                                                             * tp_value[interp.tp_pos].ars[0][j]);
 
 2840                                        Blp[4] += interp.weight
 
 2841                                                  * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 2842                                                           * tp_value[interp.tp_pos].atrs[1][j]
 
 2843                                                     + tying_point[interp.tp_pos].dN_ip[1][k]
 
 2844                                                             * tp_value[interp.tp_pos].atrs[0][j]);
 
 2847                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 2850                                          * (tp_value[interp.tp_pos].ars[0][j]
 
 2852                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2854                                                              * tying_point[interp.tp_pos]
 
 2856                                             + tp_value[interp.tp_pos].ars[1][j]
 
 2858                                                              * tying_point[interp.tp_pos]
 
 2861                                                                * tying_point[interp.tp_pos]
 
 2863                                             + tp_value[interp.tp_pos].at[j]
 
 2864                                                     * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
 
 2866                                                                * tying_point[interp.tp_pos]
 
 2870                                          * (tp_value[interp.tp_pos].atrs[0][j]
 
 2872                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2874                                                              * tying_point[interp.tp_pos]
 
 2876                                             + tp_value[interp.tp_pos].atrs[1][j]
 
 2878                                                              * tying_point[interp.tp_pos]
 
 2881                                                                * tying_point[interp.tp_pos]
 
 2886                        if (MITC_TYING_POINT::nb_ers_comp == 0) {
 
 2887                            Blp[2] = N.dN_ip[0][k] * as[j] + N.dN_ip[1][k] * ar[j];
 
 2888                            Blp[5] = N.dN_ip[0][k] * ats[j] + N.dN_ip[1][k] * atr[j];
 
 2892                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
 
 2893                                const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
 
 2895                                    if (interp.e_pos != 2) {
 
 2898                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 2899                                              * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 2902                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 2903                                              * tp_value[interp.tp_pos].atrs[interp.e_pos][j];
 
 2905                                        Blp[2] += interp.weight
 
 2906                                                  * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 2907                                                           * tp_value[interp.tp_pos].ars[1][j]
 
 2908                                                     + tying_point[interp.tp_pos].dN_ip[1][k]
 
 2909                                                             * tp_value[interp.tp_pos].ars[0][j]);
 
 2910                                        Blp[5] += interp.weight
 
 2911                                                  * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 2912                                                           * tp_value[interp.tp_pos].atrs[1][j]
 
 2913                                                     + tying_point[interp.tp_pos].dN_ip[1][k]
 
 2914                                                             * tp_value[interp.tp_pos].atrs[0][j]);
 
 2917                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 2920                                          * (tp_value[interp.tp_pos].ars[0][j]
 
 2922                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2924                                                              * tying_point[interp.tp_pos]
 
 2926                                             + tp_value[interp.tp_pos].ars[1][j]
 
 2928                                                              * tying_point[interp.tp_pos]
 
 2931                                                                * tying_point[interp.tp_pos]
 
 2933                                             + tp_value[interp.tp_pos].at[j]
 
 2934                                                     * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
 
 2936                                                                * tying_point[interp.tp_pos]
 
 2940                                          * (tp_value[interp.tp_pos].atrs[0][j]
 
 2942                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2944                                                              * tying_point[interp.tp_pos]
 
 2946                                             + tp_value[interp.tp_pos].atrs[1][j]
 
 2948                                                              * tying_point[interp.tp_pos]
 
 2951                                                                * tying_point[interp.tp_pos]
 
 2956                        if (MITC_TYING_POINT::nb_ert_comp == 0) {
 
 2957                            Blp[6] = N.dN[0][k] * at[j];
 
 2960                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
 
 2961                                const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
 
 2963                                    Blp[6] += interp.weight
 
 2964                                              * tying_point[interp.tp_pos].dN[interp.e_pos][k]
 
 2965                                              * tp_value[interp.tp_pos].at[j];
 
 2967                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 2970                                          * (tp_value[interp.tp_pos].ars[0][j]
 
 2972                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 2974                                                              * tying_point[interp.tp_pos]
 
 2976                                             + tp_value[interp.tp_pos].ars[1][j]
 
 2978                                                              * tying_point[interp.tp_pos]
 
 2981                                                                * tying_point[interp.tp_pos]
 
 2983                                             + tp_value[interp.tp_pos].at[j]
 
 2984                                                     * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
 
 2986                                                                * tying_point[interp.tp_pos]
 
 2991                        if (MITC_TYING_POINT::nb_est_comp == 0) {
 
 2992                            Blp[7] = N.dN[1][k] * at[j];
 
 2995                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
 
 2996                                const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
 
 2998                                    Blp[7] += interp.weight
 
 2999                                              * tying_point[interp.tp_pos].dN[interp.e_pos][k]
 
 3000                                              * tp_value[interp.tp_pos].at[j];
 
 3002                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3005                                          * (tp_value[interp.tp_pos].ars[0][j]
 
 3007                                                            * tying_point[interp.tp_pos].dN_ip[0][k]
 
 3009                                                              * tying_point[interp.tp_pos]
 
 3011                                             + tp_value[interp.tp_pos].ars[1][j]
 
 3013                                                              * tying_point[interp.tp_pos]
 
 3016                                                                * tying_point[interp.tp_pos]
 
 3018                                             + tp_value[interp.tp_pos].at[j]
 
 3019                                                     * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
 
 3021                                                                * tying_point[interp.tp_pos]
 
 3029                for (
int j = 0; j != 3; ++j, ++dofi1) {
 
 3030                    TP* 
const Blp = Bl[dofi1][l];
 
 3031                    Blp[0] = Blp[1] = Blp[2] = 0;
 
 3032                    if (MITC_TYING_POINT::nb_err_comp == 0) {
 
 3033                        Blp[3] = N.dN_ip[0][k] * ar[j];
 
 3036                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
 
 3037                            const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
 
 3039                                if (interp.e_pos != 2) {
 
 3040                                    Blp[3] += interp.weight
 
 3041                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 3042                                              * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3044                                    Blp[3] += interp.weight
 
 3045                                              * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 3046                                                       * tp_value[interp.tp_pos].ars[1][j]
 
 3047                                                 + tying_point[interp.tp_pos].dN_ip[1][k]
 
 3048                                                         * tp_value[interp.tp_pos].ars[0][j]);
 
 3051                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3052                                Blp[0] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3053                                          * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3054                                             + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3057                                      * (tp_value[interp.tp_pos].ars[0][j]
 
 3058                                               * (tpt[0] * tying_point[interp.tp_pos].dN_ip[0][k]
 
 3059                                                  + tpt[2] * tying_point[interp.tp_pos].dN_ip[1][k])
 
 3060                                         + tp_value[interp.tp_pos].ars[1][j]
 
 3061                                                 * (tpt[1] * tying_point[interp.tp_pos].dN_ip[1][k]
 
 3063                                                            * tying_point[interp.tp_pos]
 
 3068                    if (MITC_TYING_POINT::nb_ess_comp == 0) {
 
 3069                        Blp[4] = N.dN_ip[1][k] * as[j];
 
 3072                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
 
 3073                            const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
 
 3075                                if (interp.e_pos != 2) {
 
 3076                                    Blp[4] += interp.weight
 
 3077                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 3078                                              * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3080                                    Blp[4] += interp.weight
 
 3081                                              * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 3082                                                       * tp_value[interp.tp_pos].ars[1][j]
 
 3083                                                 + tying_point[interp.tp_pos].dN_ip[1][k]
 
 3084                                                         * tp_value[interp.tp_pos].ars[0][j]);
 
 3087                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3088                                Blp[1] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3089                                          * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3090                                             + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3093                                      * (tp_value[interp.tp_pos].ars[0][j]
 
 3094                                               * (tpt[0] * tying_point[interp.tp_pos].dN_ip[0][k]
 
 3095                                                  + tpt[2] * tying_point[interp.tp_pos].dN_ip[1][k])
 
 3096                                         + tp_value[interp.tp_pos].ars[1][j]
 
 3097                                                 * (tpt[1] * tying_point[interp.tp_pos].dN_ip[1][k]
 
 3099                                                            * tying_point[interp.tp_pos]
 
 3104                    if (MITC_TYING_POINT::nb_ers_comp == 0) {
 
 3105                        Blp[5] = N.dN_ip[0][k] * as[j] + N.dN_ip[1][k] * ar[j];
 
 3108                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
 
 3109                            const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
 
 3111                                if (interp.e_pos != 2) {
 
 3112                                    Blp[5] += interp.weight
 
 3113                                              * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
 
 3114                                              * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3116                                    Blp[5] += interp.weight
 
 3117                                              * (tying_point[interp.tp_pos].dN_ip[0][k]
 
 3118                                                       * tp_value[interp.tp_pos].ars[1][j]
 
 3119                                                 + tying_point[interp.tp_pos].dN_ip[1][k]
 
 3120                                                         * tp_value[interp.tp_pos].ars[0][j]);
 
 3123                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3124                                Blp[2] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3125                                          * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3126                                             + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3129                                      * (tp_value[interp.tp_pos].ars[0][j]
 
 3130                                               * (tpt[0] * tying_point[interp.tp_pos].dN_ip[0][k]
 
 3131                                                  + tpt[2] * tying_point[interp.tp_pos].dN_ip[1][k])
 
 3132                                         + tp_value[interp.tp_pos].ars[1][j]
 
 3133                                                 * (tpt[1] * tying_point[interp.tp_pos].dN_ip[1][k]
 
 3135                                                            * tying_point[interp.tp_pos]
 
 3140                    if (MITC_TYING_POINT::nb_ert_comp == 0) {
 
 3141                        Blp[6] = N.N[k] * ar[j];
 
 3144                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
 
 3145                            const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
 
 3147                                Blp[6] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3148                                          * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3150                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3151                                Blp[6] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3152                                          * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3153                                             + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3157                    if (MITC_TYING_POINT::nb_est_comp == 0) {
 
 3158                        Blp[7] = N.N[k] * as[j];
 
 3161                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
 
 3162                            const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
 
 3164                                Blp[7] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3165                                          * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3167                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3168                                Blp[7] += interp.weight * tying_point[interp.tp_pos].N[k]
 
 3169                                          * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3170                                             + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3175                const int j_max = nodes_type_6_dof[k] ? 3 : 2;
 
 3176                for (
int i = 3; i != 8; ++i) {
 
 3177                    const TP b0 = Bl[dofi][l][i];
 
 3178                    const TP b1 = Bl[dofi + 1][l][i];
 
 3179                    const TP b2 = Bl[dofi + 2][l][i];
 
 3180                    for (
int j = 0; j != j_max; ++j) {
 
 3181                        Bl[dofi + j][l][i] = b0 * T[k][j][0] + b1 * T[k][j][1] + b2 * T[k][j][2];
 
 3188            if (incompatible_mode) {
 
 3190                    TP* Blp = Bl[dofi][l];
 
 3191                    Blp[0] = -2 * N.IntCoor[0];
 
 3194                    std::fill(Blp + 3, Blp + 8, 0);
 
 3196                    Blp = Bl[++dofi][l];
 
 3199                    Blp[2] = -2 * N.IntCoor[1];
 
 3200                    std::fill(Blp + 3, Blp + 8, 0);
 
 3202                    Blp = Bl[++dofi][l];
 
 3205                    Blp[2] = -2 * N.IntCoor[0];
 
 3206                    std::fill(Blp + 3, Blp + 8, 0);
 
 3208                    Blp = Bl[++dofi][l];
 
 3210                    Blp[1] = -2 * N.IntCoor[1];
 
 3212                    std::fill(Blp + 3, Blp + 8, 0);
 
 3214                    TP* Blp = Bl[dofi][l];
 
 3215                    Blp[0] = 2 * N.IntCoor[0] * (-1 + 2 * N.IntCoor[0] * inc_dof[0]);
 
 3217                    Blp[2] = 4 * N.IntCoor[0] * N.IntCoor[1] * inc_dof[1];
 
 3218                    std::fill(Blp + 3, Blp + 8, 0);
 
 3220                    Blp = Bl[++dofi][l];
 
 3222                    Blp[1] = 4 * N.IntCoor[1] * N.IntCoor[1] * inc_dof[1];
 
 3223                    Blp[2] = -2 * N.IntCoor[1] + 4 * N.IntCoor[0] * N.IntCoor[1] * inc_dof[0];
 
 3224                    std::fill(Blp + 3, Blp + 8, 0);
 
 3226                    Blp = Bl[++dofi][l];
 
 3227                    Blp[0] = 4 * N.IntCoor[0] * N.IntCoor[0] * inc_dof[2];
 
 3229                    Blp[2] = -2 * N.IntCoor[0] + 4 * N.IntCoor[0] * N.IntCoor[1] * inc_dof[3];
 
 3230                    std::fill(Blp + 3, Blp + 8, 0);
 
 3232                    Blp = Bl[++dofi][l];
 
 3234                    Blp[1] = 2 * N.IntCoor[1] * (-1 + 2 * N.IntCoor[1] * inc_dof[3]);
 
 3235                    Blp[2] = 4 * N.IntCoor[0] * N.IntCoor[1] * inc_dof[2];
 
 3236                    std::fill(Blp + 3, Blp + 8, 0);
 
 3242        if ((d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]) && !linear) {
 
 3243            for (
int k = 0; k != nb_nodes; ++k) {
 
 3244                for (
int j = 0; j != 3; ++j) {
 
 3245                    if (MITC_TYING_POINT::nb_err_comp == 0) {
 
 3246                        h[k][j] += S[l][3] * N.dN[0][k] * ar[j];
 
 3248                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
 
 3249                            const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
 
 3251                                if (interp.e_pos != 2) {
 
 3252                                    h[k][j] += S[l][3] * interp.weight
 
 3253                                               * tying_point[interp.tp_pos].dN[interp.e_pos][k]
 
 3254                                               * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3256                                    h[k][j] += S[l][3] * interp.weight
 
 3257                                               * (tying_point[interp.tp_pos].dN[0][k]
 
 3258                                                        * tp_value[interp.tp_pos].ars[1][j]
 
 3259                                                  + tying_point[interp.tp_pos].dN[1][k]
 
 3260                                                          * tp_value[interp.tp_pos].ars[0][j]);
 
 3263                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3265                                      S[l][3] * interp.weight
 
 3266                                      * (tp_value[interp.tp_pos].ars[0][j]
 
 3267                                               * (tpt[0] * tying_point[interp.tp_pos].dN[0][k]
 
 3268                                                  + tpt[2] * tying_point[interp.tp_pos].dN[1][k])
 
 3269                                         + tp_value[interp.tp_pos].ars[1][j]
 
 3270                                                 * (tpt[1] * tying_point[interp.tp_pos].dN[1][k]
 
 3272                                                            * tying_point[interp.tp_pos].dN[0][k]));
 
 3276                    if (MITC_TYING_POINT::nb_ess_comp == 0) {
 
 3277                        h[k][j] += S[l][4] * N.dN[1][k] * as[j];
 
 3279                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
 
 3280                            const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
 
 3282                                if (interp.e_pos != 2) {
 
 3283                                    h[k][j] += S[l][4] * interp.weight
 
 3284                                               * tying_point[interp.tp_pos].dN[interp.e_pos][k]
 
 3285                                               * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3287                                    h[k][j] += S[l][4] * interp.weight
 
 3288                                               * (tying_point[interp.tp_pos].dN[0][k]
 
 3289                                                        * tp_value[interp.tp_pos].ars[1][j]
 
 3290                                                  + tying_point[interp.tp_pos].dN[1][k]
 
 3291                                                          * tp_value[interp.tp_pos].ars[0][j]);
 
 3294                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3296                                      S[l][4] * interp.weight
 
 3297                                      * (tp_value[interp.tp_pos].ars[0][j]
 
 3298                                               * (tpt[0] * tying_point[interp.tp_pos].dN[0][k]
 
 3299                                                  + tpt[2] * tying_point[interp.tp_pos].dN[1][k])
 
 3300                                         + tp_value[interp.tp_pos].ars[1][j]
 
 3301                                                 * (tpt[1] * tying_point[interp.tp_pos].dN[1][k]
 
 3303                                                            * tying_point[interp.tp_pos].dN[0][k]));
 
 3307                    if (MITC_TYING_POINT::nb_ers_comp == 0) {
 
 3308                        h[k][j] += S[l][5] * (N.dN[0][k] * as[j] + N.dN[1][k] * ar[j]);
 
 3310                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
 
 3311                            const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
 
 3313                                if (interp.e_pos != 2) {
 
 3314                                    h[k][j] += S[l][5] * interp.weight
 
 3315                                               * tying_point[interp.tp_pos].dN[interp.e_pos][k]
 
 3316                                               * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3318                                    h[k][j] += S[l][5] * interp.weight
 
 3319                                               * (tying_point[interp.tp_pos].dN[0][k]
 
 3320                                                        * tp_value[interp.tp_pos].ars[1][j]
 
 3321                                                  + tying_point[interp.tp_pos].dN[1][k]
 
 3322                                                          * tp_value[interp.tp_pos].ars[0][j]);
 
 3325                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3327                                      S[l][5] * interp.weight
 
 3328                                      * (tp_value[interp.tp_pos].ars[0][j]
 
 3329                                               * (tpt[0] * tying_point[interp.tp_pos].dN[0][k]
 
 3330                                                  + tpt[2] * tying_point[interp.tp_pos].dN[1][k])
 
 3331                                         + tp_value[interp.tp_pos].ars[1][j]
 
 3332                                                 * (tpt[1] * tying_point[interp.tp_pos].dN[1][k]
 
 3334                                                            * tying_point[interp.tp_pos].dN[0][k]));
 
 3338                    if (MITC_TYING_POINT::nb_ert_comp == 0) {
 
 3339                        h[k][j] += S[l][6] * N.N[k] * ar[j];
 
 3341                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
 
 3342                            const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
 
 3344                                h[k][j] += S[l][6] * interp.weight * tying_point[interp.tp_pos].N[k]
 
 3345                                           * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3347                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3348                                h[k][j] += S[l][6] * interp.weight * tying_point[interp.tp_pos].N[k]
 
 3349                                           * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3350                                              + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3354                    if (MITC_TYING_POINT::nb_est_comp == 0) {
 
 3355                        h[k][j] += S[l][7] * N.N[k] * as[j];
 
 3357                        for (
int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
 
 3358                            const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
 
 3360                                h[k][j] += S[l][7] * interp.weight * tying_point[interp.tp_pos].N[k]
 
 3361                                           * tp_value[interp.tp_pos].ars[interp.e_pos][j];
 
 3363                                const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3364                                h[k][j] += S[l][7] * interp.weight * tying_point[interp.tp_pos].N[k]
 
 3365                                           * (tpt[3] * tp_value[interp.tp_pos].ars[0][j]
 
 3366                                              + tpt[4] * tp_value[interp.tp_pos].ars[1][j]);
 
 3374        if (!d_value_d_dofdotdot_null || compute_value_from_dofdotdot) {
 
 3378            const int number_of_layer = 
property->number_of_layer();
 
 3380                  b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin,
 
 3392                TP inertia_force[3];
 
 3394            std::vector<Data> data_vec;
 
 3397            for (
int layer = 0; layer != number_of_layer; ++layer) {
 
 3401                      b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
 
 3403                for (
int lt = 0; lt != nb_tt_integration; ++lt) {
 
 3404                    data_vec.push_back(Data());
 
 3405                    Data& 
data = data_vec.back();
 
 3406                    data.omega = l_begin + 0.5 * (1 + n_omega[lt]) * (l_end - l_begin);
 
 3407                    for (
int j = 0; j != 3; ++j) {
 
 3408                        data.G[0][j] = Ar[j] + 
data.omega * Atr[j];
 
 3409                        data.G[1][j] = As[j] + 
data.omega * Ats[j];
 
 3410                        data.G[2][j] = At[j];
 
 3412                    data.density = 
property->get_density(layer) * gauss_point[l].weight
 
 3414                    std::fill_n(
data.inertia_force, 3, 0.);
 
 3417            if (non_structural_mass != 0) {
 
 3418                data_vec.push_back(Data());
 
 3419                Data& 
data = data_vec.back();
 
 3421                for (
int j = 0; j != 3; ++j) {
 
 3422                    data.G[0][j] = Ar[j] + 
data.omega * Atr[j];
 
 3423                    data.G[1][j] = As[j] + 
data.omega * Ats[j];
 
 3424                    data.G[2][j] = At[j];
 
 3428                std::fill_n(
data.inertia_force, 3, 0.);
 
 3434            for (
auto data : data_vec) {
 
 3435                TP Nl[max_nb_dof][3];
 
 3436                std::fill_n(Nl[0], number_of_dof * 3, 0);
 
 3438                for (
int k = 0; k != nb_nodes; ++k) {
 
 3439                    const TP Nlk = N.N[k];
 
 3440                    if (SHAPE_IP_OOP_SAME) {
 
 3441                        Nl[dofi][0] = Nl[dofi + 1][1] = Nl[dofi + 2][2] = Nlk;
 
 3443                    } 
else if (k < nb_nodes_ip) {
 
 3444                        Nl[dofi][0] = Nl[dofi + 1][1] = Nl[dofi + 2][2] = N.N_ip[k];
 
 3447                    const TP Nlk_omega = Nlk * 
data.omega;
 
 3448                    if (nodes_type_6_dof[k]) {
 
 3450                        Nl[dofi][1] = Nlk_omega * -D[k][2];
 
 3451                        Nl[dofi][2] = Nlk_omega * D[k][1];
 
 3453                        Nl[++dofi][0] = Nlk_omega * D[k][2];
 
 3455                        Nl[dofi][2] = Nlk_omega * -D[k][0];
 
 3457                        Nl[++dofi][0] = Nlk_omega * -D[k][1];
 
 3458                        Nl[dofi][1] = Nlk_omega * D[k][0];
 
 3462                            const TP* 
const Dk = Dr[k][0];
 
 3463                            Nl[dofi][0] = Nlk_omega * (D[k][2] * Dk[1] - D[k][1] * Dk[2]);
 
 3464                            Nl[dofi][1] = Nlk_omega * (-D[k][2] * Dk[0] + D[k][0] * Dk[2]);
 
 3465                            Nl[dofi][2] = Nlk_omega * (D[k][1] * Dk[0] - D[k][0] * Dk[1]);
 
 3468                            const TP* 
const Dk = Dr[k][1];
 
 3469                            Nl[++dofi][0] = Nlk_omega * (D[k][2] * Dk[1] - D[k][1] * Dk[2]);
 
 3470                            Nl[dofi][1] = Nlk_omega * (-D[k][2] * Dk[0] + D[k][0] * Dk[2]);
 
 3471                            Nl[dofi][2] = Nlk_omega * (D[k][1] * Dk[0] - D[k][0] * Dk[1]);
 
 3476                if (incompatible_mode) {
 
 3478                    const double r = 1 - N.IntCoor[0] * N.IntCoor[0];
 
 3479                    const double s = 1 - N.IntCoor[1] * N.IntCoor[1];
 
 3480                    for (
int j = 0; j != 3; ++j) {
 
 3481                        Nl[dofi][j] = r * 
data.G[0][j];
 
 3482                        Nl[dofi + 1][j] = s * 
data.G[0][j];
 
 3483                        Nl[dofi + 2][j] = r * 
data.G[1][j];
 
 3484                        Nl[dofi + 3][j] = s * 
data.G[1][j];
 
 3487                if (compute_value_from_dofdotdot) {
 
 3489                          'N', 3, number_of_dof, 
data.density, Nl[0], 3, &dof(0, 2), 1, 1.0,
 
 3490                          data.inertia_force, 1);
 
 3492                          'T', 3, number_of_dof, 1.0, Nl[0], 3, 
data.inertia_force, 1, 1.0,
 
 3495                if (!d_value_d_dofdotdot_null) {
 
 3496                    TP* ptr = &d_value_d_dof[2](0, 0);
 
 3497                    for (
int j = 0; j != 3; ++j) {
 
 3498                        blas::spr(
'U', number_of_dof, 
data.density, Nl[0] + j, 3, ptr);
 
 3509    const TP drill_stiffness_factor = 
property->get_drill_stiffness();
 
 3510    TP drill_stiffness_gauss[nb_gauss] = {};
 
 3511    if (drill_stiffness_factor != 0) {
 
 3512        for (
int l = 0; l != nb_gauss; ++l) {
 
 3513            const OnePoint& N = gauss_point[l];
 
 3517                  b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin,
 
 3524            for (
int k = 0; k != nb_nodes; ++k) {
 
 3525                for (
int j = 0; j != 3; ++j) {
 
 3526                    Ar[j] += N.dN[0][k] * R[k][j];
 
 3527                    As[j] += N.dN[1][k] * R[k][j];
 
 3536                area = 
norm_3(normal) * gauss_point[l].weight;
 
 3539            const int number_of_layer = 
property->number_of_layer();
 
 3540            for (
int layer = 0; layer != number_of_layer; ++layer) {
 
 3544                      b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
 
 3546                drill_stiffness_gauss[l] += 
property->get_inplane_shear_modulus(layer) * area
 
 3547                                            * (l_end - l_begin) * drill_stiffness_factor;
 
 3554    if (!value.is_null()) {
 
 3556              'T', 8 * nb_gauss, number_of_dof, 1.0, Bl[0][0], 8 * nb_gauss, S[0], 1,
 
 3557              (compute_value_from_dofdotdot ? 1.0 : 0.0), &value[0], 1);
 
 3558        if (drill_stiffness_factor != 0) {
 
 3560            for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 3562                if (nodes_type_6_dof[k]) {
 
 3563                    if (drill_stabilized_node[k]) {
 
 3564                        for (
int l = 0; l != nb_gauss; ++l) {
 
 3565                            const OnePointShape& N = gauss_point[l];
 
 3566                            const TP drill = drill_stiffness_gauss[l] * 
dot_3(w_dof[k], D[k])
 
 3568                            for (
int i = 0; i != 3; ++i) { value[p + i] += drill * D[k][i]; }
 
 3580    if (d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]) {
 
 3585            TP* out_ptr = &d_value_d_dof[0](0, 0);
 
 3586            for (
int i = 0; i != number_of_dof;) {
 
 3587                TP tmp[nb_gauss][8] = {};
 
 3588                for (
int l = 0; l != nb_gauss; ++l) {
 
 3590                        const TP* cc = C[l];
 
 3591                        for (
int jj = 0; jj != 8; ++jj) {
 
 3592                            tmp[l][jj] += *cc * Bl[i][l][jj];
 
 3594                            for (
int ii = jj + 1; ii < 8; ++ii, ++cc) {
 
 3595                                tmp[l][jj] += *cc * Bl[i][l][ii];
 
 3596                                tmp[l][ii] += *cc * Bl[i][l][jj];
 
 3603                      'T', nb_gauss * 8, i, 1.0, Bl[0][0], nb_gauss * 8, tmp[0], 1, 0.0, out_ptr,
 
 3612            for (
int kj = 0; kj != nb_nodes; ++kj) {
 
 3613                const int j_max = nodes_type_6_dof[kj] ? 3 : 2;
 
 3615                for (
int ki = 0; ki != nb_nodes; ++ki) {
 
 3616                    const int i_max = nodes_type_6_dof[ki] ? 3 : 2;
 
 3620                    for (
int l = 0; l != nb_gauss; ++l) {
 
 3621                        const OnePoint& N = gauss_point[l];
 
 3623                        if (MITC_TYING_POINT::nb_err_comp == 0) {
 
 3624                            const double tmp = N.dN_ip[0][ki] * N.dN_ip[0][kj];
 
 3629                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
 
 3630                                const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
 
 3632                                    if (interp.e_pos != 2) {
 
 3633                                        tmp += interp.weight
 
 3634                                               * tying_point[interp.tp_pos].dN_ip[interp.e_pos][ki]
 
 3635                                               * tying_point[interp.tp_pos].dN_ip[interp.e_pos][kj];
 
 3639                                              * (tying_point[interp.tp_pos].dN_ip[0][ki]
 
 3640                                                       * tying_point[interp.tp_pos].dN_ip[1][kj]
 
 3641                                                 + tying_point[interp.tp_pos].dN_ip[1][ki]
 
 3642                                                         * tying_point[interp.tp_pos].dN_ip[0][kj]);
 
 3645                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3646                                    tmp += interp.weight
 
 3647                                           * (tying_point[interp.tp_pos].dN_ip[0][ki]
 
 3649                                                             * tying_point[interp.tp_pos]
 
 3652                                                               * tying_point[interp.tp_pos]
 
 3654                                              + tying_point[interp.tp_pos].dN_ip[0][kj]
 
 3656                                                               * tying_point[interp.tp_pos]
 
 3659                                                                 * tying_point[interp.tp_pos]
 
 3666                        if (MITC_TYING_POINT::nb_ess_comp == 0) {
 
 3667                            const double tmp = N.dN_ip[1][ki] * N.dN_ip[1][kj];
 
 3672                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
 
 3673                                const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
 
 3675                                    if (interp.e_pos != 2) {
 
 3676                                        tmp += interp.weight
 
 3677                                               * tying_point[interp.tp_pos].dN_ip[interp.e_pos][ki]
 
 3678                                               * tying_point[interp.tp_pos].dN_ip[interp.e_pos][kj];
 
 3682                                              * (tying_point[interp.tp_pos].dN_ip[0][ki]
 
 3683                                                       * tying_point[interp.tp_pos].dN_ip[1][kj]
 
 3684                                                 + tying_point[interp.tp_pos].dN_ip[1][ki]
 
 3685                                                         * tying_point[interp.tp_pos].dN_ip[0][kj]);
 
 3688                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3689                                    tmp += interp.weight
 
 3690                                           * (tying_point[interp.tp_pos].dN_ip[0][ki]
 
 3692                                                             * tying_point[interp.tp_pos]
 
 3695                                                               * tying_point[interp.tp_pos]
 
 3697                                              + tying_point[interp.tp_pos].dN_ip[0][kj]
 
 3699                                                               * tying_point[interp.tp_pos]
 
 3702                                                                 * tying_point[interp.tp_pos]
 
 3709                        if (MITC_TYING_POINT::nb_ers_comp == 0) {
 
 3711                                  N.dN_ip[0][ki] * N.dN_ip[1][kj] + N.dN_ip[1][ki] * N.dN_ip[0][kj];
 
 3716                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
 
 3717                                const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
 
 3719                                    if (interp.e_pos != 2) {
 
 3720                                        tmp += interp.weight
 
 3721                                               * tying_point[interp.tp_pos].dN_ip[interp.e_pos][ki]
 
 3722                                               * tying_point[interp.tp_pos].dN_ip[interp.e_pos][kj];
 
 3726                                              * (tying_point[interp.tp_pos].dN_ip[0][ki]
 
 3727                                                       * tying_point[interp.tp_pos].dN_ip[1][kj]
 
 3728                                                 + tying_point[interp.tp_pos].dN_ip[1][ki]
 
 3729                                                         * tying_point[interp.tp_pos].dN_ip[0][kj]);
 
 3732                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
 
 3733                                    tmp += interp.weight
 
 3734                                           * (tying_point[interp.tp_pos].dN_ip[0][ki]
 
 3736                                                             * tying_point[interp.tp_pos]
 
 3739                                                               * tying_point[interp.tp_pos]
 
 3741                                              + tying_point[interp.tp_pos].dN_ip[0][kj]
 
 3743                                                               * tying_point[interp.tp_pos]
 
 3746                                                                 * tying_point[interp.tp_pos]
 
 3755                        if (MITC_TYING_POINT::nb_ert_comp == 0) {
 
 3756                            qrs = S[l][6] * N.dN_ip[0][ki] * N.N_ip[kj];
 
 3757                            qsr = S[l][6] * N.dN_ip[0][kj] * N.N_ip[ki];
 
 3759                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
 
 3760                                const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
 
 3762                                    qrs += S[l][6] * interp.weight
 
 3763                                           * tying_point[interp.tp_pos].dN[interp.e_pos][ki]
 
 3764                                           * tying_point[interp.tp_pos].N[kj];
 
 3765                                    qsr += S[l][6] * interp.weight
 
 3766                                           * tying_point[interp.tp_pos].dN[interp.e_pos][kj]
 
 3767                                           * tying_point[interp.tp_pos].N[ki];
 
 3769                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3770                                    qrs += S[l][6] * interp.weight
 
 3771                                           * tying_point[interp.tp_pos].N[kj]
 
 3772                                           * (tpt[3] * tying_point[interp.tp_pos].dN[0][ki]
 
 3773                                              + tpt[4] * tying_point[interp.tp_pos].dN[1][ki]);
 
 3774                                    qsr += S[l][6] * interp.weight
 
 3775                                           * tying_point[interp.tp_pos].N[ki]
 
 3776                                           * (tpt[3] * tying_point[interp.tp_pos].dN[0][kj]
 
 3777                                              + tpt[4] * tying_point[interp.tp_pos].dN[1][kj]);
 
 3781                        if (MITC_TYING_POINT::nb_est_comp == 0) {
 
 3782                            qrs += S[l][7] * N.dN_ip[1][ki] * N.N_ip[kj];
 
 3783                            qsr += S[l][7] * N.dN_ip[1][kj] * N.N_ip[ki];
 
 3785                            for (
int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
 
 3786                                const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
 
 3788                                    qrs += S[l][7] * interp.weight
 
 3789                                           * tying_point[interp.tp_pos].dN[interp.e_pos][ki]
 
 3790                                           * tying_point[interp.tp_pos].N[kj];
 
 3791                                    qsr += S[l][7] * interp.weight
 
 3792                                           * tying_point[interp.tp_pos].dN[interp.e_pos][kj]
 
 3793                                           * tying_point[interp.tp_pos].N[ki];
 
 3795                                    const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
 
 3796                                    qrs += S[l][7] * interp.weight
 
 3797                                           * tying_point[interp.tp_pos].N[kj]
 
 3798                                           * (tpt[3] * tying_point[interp.tp_pos].dN[0][ki]
 
 3799                                              + tpt[4] * tying_point[interp.tp_pos].dN[1][ki]);
 
 3800                                    qsr += S[l][7] * interp.weight
 
 3801                                           * tying_point[interp.tp_pos].N[ki]
 
 3802                                           * (tpt[3] * tying_point[interp.tp_pos].dN[0][kj]
 
 3803                                              + tpt[4] * tying_point[interp.tp_pos].dN[1][kj]);
 
 3812                    for (
int i = 0; i != i_max; ++i) {
 
 3813                        for (
int j = 0; j != 3; ++j) {
 
 3814                            d_value_d_dof[0](VD_row + i + 3, VD_col + j) += m_qsr * T[ki][i][j];
 
 3817                    for (
int i = 0; i != 3; ++i) {
 
 3818                        for (
int j = 0; j != j_max; ++j) {
 
 3819                            d_value_d_dof[0](VD_row + i, VD_col + j + 3) += m_qrs * T[kj][j][i];
 
 3823                        for (
int i = 0; i != 3; ++i) {
 
 3824                            d_value_d_dof[0](VD_row + i, VD_col + i) += n;
 
 3831                              bki[0] * w_dof[ki][0] + bki[1] * w_dof[ki][1] + bki[2] * w_dof[ki][2];
 
 3835                        if (norm_w[ki] < 0.001) {
 
 3836                            const TP norm_w2 = norm_w[ki] * norm_w[ki];
 
 3837                            c10 = 1.0 / 6 + norm_w2 / 180;
 
 3838                            c3 = 1.0 / 6 + norm_w2 / 360;
 
 3839                            c11 = -1.0 / 360 - norm_w2 / 7560;
 
 3841                            TP sin_w = std::sin(norm_w[ki]);
 
 3842                            TP cos_w_1 = std::sin(norm_w[ki] / 2);
 
 3843                            cos_w_1 = -2 * cos_w_1 * cos_w_1;
 
 3844                            c10 = (sin_w - norm_w[ki]) / (2 * norm_w[ki] * cos_w_1);
 
 3845                            TP norm_w2 = norm_w[ki] * norm_w[ki];
 
 3846                            c11 = (4 * cos_w_1 + norm_w2 + norm_w[ki] * sin_w)
 
 3847                                  / (2 * norm_w2 * norm_w2 * cos_w_1);
 
 3848                            c3 = (norm_w[ki] * sin_w + 2 * cos_w_1) / (norm_w2 * cos_w_1);
 
 3852                              - (d[ki][0] * h[ki][0] + d[ki][1] * h[ki][1] + d[ki][2] * h[ki][2]);
 
 3854                        for (
int i = 0; i != 3; ++i) { bki[i] = -c3 * bki[i] + c11 * w_dof[ki][i]; }
 
 3856                        for (
int i = 0; i != 3; ++i) {
 
 3857                            for (
int j = 0; j != 3; ++j) {
 
 3859                                          * (d[ki][i] * h[ki][j] + d[ki][j] * h[ki][i]
 
 3860                                             + bki[i] * w_dof[ki][j] + bki[j] * w_dof[ki][i]);
 
 3864                        for (
int j = 0; j != i_max; ++j) {
 
 3866                            for (
int i = 0; i != 3; ++i) {
 
 3867                                MHT3[i] = M[0][i] * HT3[ki][j][0] + M[1][i] * HT3[ki][j][1]
 
 3868                                          + M[2][i] * HT3[ki][j][2];
 
 3870                            for (
int i = j; i != i_max; ++i) {
 
 3871                                d_value_d_dof[0](VD_row + i + 3, VD_col + j + 3) +=
 
 3872                                      HT3[ki][i][0] * MHT3[0] + HT3[ki][i][1] * MHT3[1]
 
 3873                                      + HT3[ki][i][2] * MHT3[2];
 
 3877                    VD_row += 3 + i_max;
 
 3879                VD_col += 3 + j_max;
 
 3881            if (incompatible_mode) {
 
 3889                for (
int l = 0; l != nb_gauss; ++l) {
 
 3890                    const OnePoint& N = gauss_point[l];
 
 3891                    const double nrr = N.IntCoor[0] * N.IntCoor[0];
 
 3892                    const double nss = N.IntCoor[1] * N.IntCoor[1];
 
 3893                    const double nrs = N.IntCoor[0] * N.IntCoor[1];
 
 3894                    a00 += nrr * S[l][0];
 
 3895                    a01 += nrs * S[l][2];
 
 3896                    a11 += nss * S[l][1];
 
 3897                    a22 += nrr * S[l][0];
 
 3898                    a23 += nrs * S[l][2];
 
 3899                    a33 += nss * S[l][1];
 
 3901                d_value_d_dof[0](VD_col, VD_col) += 4 * a00;
 
 3902                d_value_d_dof[0](VD_col, VD_col + 1) += 4 * a01;
 
 3903                d_value_d_dof[0](VD_col + 1, VD_col + 1) += 4 * a11;
 
 3904                d_value_d_dof[0](VD_col + 2, VD_col + 2) += 4 * a22;
 
 3905                d_value_d_dof[0](VD_col + 2, VD_col + 3) += 4 * a23;
 
 3906                d_value_d_dof[0](VD_col + 3, VD_col + 3) += 4 * a33;
 
 3910        if (drill_stiffness_factor != 0) {
 
 3912            for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 3914                if (nodes_type_6_dof[k]) {
 
 3915                    if (drill_stabilized_node[k]) {
 
 3916                        for (
int l = 0; l != nb_gauss; ++l) {
 
 3917                            const OnePointShape& N = gauss_point[l];
 
 3918                            for (
int i = 0; i != 3; ++i) {
 
 3919                                for (
int j = i; j != 3; ++j) {
 
 3920                                    d_value_d_dof[0](p + i, p + j) += drill_stiffness_gauss[l]
 
 3921                                                                      * D[k][i] * D[k][j] * N.N[k]
 
 3937      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 3938      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 3939      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 3940void b2000::ShellMITC<
 
 3941      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 3942      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::
 
 3943      field_edge_integration(
 
 3944            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dof,
 
 3945            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdot,
 
 3946            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdotdot,
 
 3947            const Field<TP>& field, 
int edge, b2linalg::Index& dof_numbering,
 
 3948            b2linalg::Vector<TP, b2linalg::Vdense>& discretised_field,
 
 3949            b2linalg::Matrix<TP, b2linalg::Mpacked>& d_discretised_field_d_dof) {
 
 3954    if (edge < 0 || edge >= SHAPE::num_edge) {
 
 3955        Exception() << 
"Invalid edge number " << edge + 1 << 
" for element " 
 3956                    << this->get_object_name() << 
"." << 
THROW;
 
 3961    if (discretised_field.is_null()) { 
return; }
 
 3963    typedef typename SHAPE::edge_shape_type_0 edge_shape;
 
 3967    size_t edn[3 * edge_shape::num_nodes];
 
 3968    for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 3969        const size_t ii = SHAPE::edge_node[edge][i];
 
 3970        assert(ii < (
size_t)nb_nodes_ip);
 
 3972        if (nodes_type_6_dof[ii]) {
 
 3973            dof6_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
 
 3975            dof5_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
 
 3977        std::copy(ndn, ndn + 3, edn + 3 * i);
 
 3982    b2linalg::Index adn;
 
 3983    b2linalg::Index dpos;
 
 3984    for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 3985        const size_t ii = SHAPE::edge_node[edge][i];
 
 3986        assert(ii < (
size_t)nb_nodes_ip);
 
 3988        dpos.push_back(adn.size());
 
 3989        if (nodes_type_6_dof[ii]) {
 
 3990            dof6_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
 
 3991            adn.insert(adn.end(), ndn, ndn + 6);
 
 3993            dof5_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
 
 3994            adn.insert(adn.end(), ndn, ndn + 5);
 
 3997    dpos.push_back(adn.size());
 
 3998    assert(dpos.size() == edge_shape::num_nodes + 1);
 
 3999    assert(dpos.front() == 0 && dpos.back() == adn.size());
 
 4001    discretised_field.resize(adn.size());
 
 4002    discretised_field.set_zero();
 
 4003    if (!dof_numbering.is_null()) { dof_numbering = adn; }
 
 4007    double enode_xi[edge_shape::num_nodes][3] = {};
 
 4008    TP enode_x_0[edge_shape::num_nodes + 1][3] = {};
 
 4009    TP enode_x_d[edge_shape::num_nodes][3] = {};
 
 4010    for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 4011        const size_t ii = SHAPE::edge_node[edge][i];
 
 4012        enode_xi[i][0] = SHAPE::r_node[ii];
 
 4013        enode_xi[i][1] = SHAPE::s_node[ii];
 
 4014        enode_xi[i][2] = 0.0;
 
 4015        if (nodes_type_6_dof[ii]) {
 
 4016            dof6_node_type::get_coor_s(enode_x_0[i], nodes[ii]);
 
 4018            dof5_node_type::get_coor_s(enode_x_0[i], nodes[ii]);
 
 4020        std::copy(enode_x_0[i], enode_x_0[i + 1], enode_x_d[i]);
 
 4021        if (!dof.is_null()) {
 
 4022            for (
int j = 0; j != 3; ++j) { enode_x_d[i][j] += dof[edn[3 * i + j]]; }
 
 4028    ischeme.set_order(edge_shape::order);
 
 4029    for (
int l = 0; l != ischeme.get_num(); ++l) {
 
 4032        ischeme.get_point(l, r, weight);
 
 4034        double N[edge_shape::num_nodes];
 
 4035        edge_shape::eval_h(&r, N);
 
 4036        double dN[edge_shape::num_nodes][1];
 
 4037        edge_shape::eval_h_derived(&r, dN);
 
 4042        for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 4043            for (
int j = 0; j != 3; ++j) {
 
 4044                dr_0[j] += dN[i][0] * enode_x_0[i][j];
 
 4045                dr_d[j] += dN[i][0] * enode_x_d[i][j];
 
 4058        for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 4059            for (
int j = 0; j != 3; ++j) {
 
 4060                xi[j] += N[i] * enode_xi[i][j];
 
 4061                x_0[j] += N[i] * enode_x_0[i][j];
 
 4062                x_d[j] += N[i] * enode_x_d[i][j];
 
 4067        b2linalg::Vector<TP, b2linalg::Vdense> value;
 
 4069              b2linalg::Vector<double, b2linalg::Vdense_constref>(3, xi),
 
 4070              b2linalg::Vector<TP, b2linalg::Vdense_constref>(3, x_0),
 
 4071              b2linalg::Vector<TP, b2linalg::Vdense_constref>(3, x_d),
 
 4072              b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>(3, 1, dr_0),
 
 4073              b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>(3, 1, dr_d), value);
 
 4077        if (line_load_as_moment) {
 
 4078            for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 4079                for (
size_t j = 3; j < dpos[i + 1] - dpos[i]; ++j) {
 
 4080                    discretised_field[dpos[i] + j] += length * value[j - 3] * N[i];
 
 4084            for (
int i = 0; i != edge_shape::num_nodes; ++i) {
 
 4085                for (
size_t j = 0; j != 3; ++j) {
 
 4086                    discretised_field[dpos[i] + j] += length * value[j] * N[i];
 
 4094      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 4095      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 4096      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 4097void b2000::ShellMITC<
 
 4098      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 4099      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::
 
 4100      field_face_integration(
 
 4101            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dof,
 
 4102            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdot,
 
 4103            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdotdot,
 
 4104            const Field<TP>& field, 
int face, b2linalg::Index& dof_numbering,
 
 4105            b2linalg::Vector<TP, b2linalg::Vdense>& discretised_field,
 
 4106            b2linalg::Matrix<TP, b2linalg::Mpacked>& d_discretised_field_d_dof) {
 
 4115                             << 
" cannot integrate face/edge fields over edges" << 
THROW;
 
 4121    TP D[nb_nodes][3] = {};
 
 4122    TP Dd[nb_nodes][3] = {};
 
 4125    TP Dr[nb_nodes][2][3] = {};
 
 4128    TP R[nb_nodes][6] = {};
 
 4130    int number_of_dof = incompatible_mode ? 4 : 0;
 
 4132    for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 4133        if (nodes_type_6_dof[k]) {
 
 4135            dof6_node_type::get_coor_s(R[k], nodes[k]);
 
 4138            dof5_node_type::get_coor_s(R[k], nodes[k]);
 
 4142    for (
int k = nb_nodes_ip; k != nb_nodes; ++k) {
 
 4144        dof2_node_type::get_coor_s(R[k], nodes[k]);
 
 4148    DIRECTOR_ESTIMATION::compute_normal(R, D);
 
 4149    for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 4150        if (!nodes_type_6_dof[k]
 
 4151            && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
 4152            if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
 4155            std::copy(R[k] + 3, R[k] + 6, D[k]);
 
 4160    for (
int k = 0; k != nb_nodes; ++k) {
 
 4161        if (!nodes_type_6_dof[k]) {
 
 4162            const TP e2[3] = {0, 1, 0};
 
 4165                const TP e2[3] = {1, 0, 0};
 
 4175    TP u_dof[nb_nodes][3];
 
 4178    TP w_dof[nb_nodes][3];
 
 4184    if (dof.is_null()) {
 
 4185        std::fill_n(u_dof[0], nb_nodes * 3, 0);
 
 4186        std::fill_n(w_dof[0], nb_nodes * 3, 0);
 
 4187        if (incompatible_mode) { std::fill_n(inc_dof, 4, 0); }
 
 4189        const TP* alpha = &dof[0];
 
 4190        for (
int k = 0; k != nb_nodes; ++k) {
 
 4191            if (nodes_type_6_dof[k]) {
 
 4192                std::copy(alpha, alpha + 3, u_dof[k]);
 
 4194                std::copy(alpha, alpha + 3, w_dof[k]);
 
 4197                std::copy(alpha, alpha + 3, u_dof[k]);
 
 4198                w_dof[k][0] = Dr[k][0][0] * alpha[3] + Dr[k][1][0] * alpha[4];
 
 4199                w_dof[k][1] = Dr[k][0][1] * alpha[3] + Dr[k][1][1] * alpha[4];
 
 4200                w_dof[k][2] = Dr[k][0][2] * alpha[3] + Dr[k][1][2] * alpha[4];
 
 4211    TP norm_w[nb_nodes];
 
 4215    for (
int k = 0; k != nb_nodes; ++k) {
 
 4216        const TP* w = w_dof[k];
 
 4217        const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
 4218        norm_w[k] = std::sqrt(ww);
 
 4221            const TP w0_2 = w[0] * w[0];
 
 4222            const TP w1_2 = w[1] * w[1];
 
 4223            const TP w2_2 = w[2] * w[2];
 
 4224            O2[0] = -w2_2 - w1_2;
 
 4225            O2[1] = w[0] * w[1];
 
 4226            O2[2] = w[0] * w[2];
 
 4227            O2[3] = -w2_2 - w0_2;
 
 4228            O2[4] = w[1] * w[2];
 
 4229            O2[5] = -w1_2 - w0_2;
 
 4238            s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
 4239            c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
 4241            s = std::sin(norm_w[k]) / norm_w[k];
 
 4242            c = std::sin(norm_w[k] * 0.5) / norm_w[k];
 
 4246            TP* 
const dk = d[k];
 
 4247            const TP* 
const Dk = D[k];
 
 4249                TP* 
const Ddk = Dd[k];
 
 4250                Ddk[0] = Dk[2] * w_dof[k][1] - Dk[1] * w_dof[k][2];
 
 4251                Ddk[1] = -Dk[2] * w_dof[k][0] + Dk[0] * w_dof[k][2];
 
 4252                Ddk[2] = Dk[1] * w_dof[k][0] - Dk[0] * w_dof[k][1];
 
 4253                dk[0] = Ddk[0] + Dk[0];
 
 4254                dk[1] = Ddk[1] + Dk[1];
 
 4255                dk[2] = Ddk[2] + Dk[2];
 
 4257                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 4258                        + (s * w[1] + c * O2[2]) * Dk[2];
 
 4259                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 4260                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
 4261                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 4262                        + (1 + c * O2[5]) * Dk[2];
 
 4263                TP* 
const Ddk = Dd[k];
 
 4264                Ddk[0] = dk[0] - Dk[0];
 
 4265                Ddk[1] = dk[1] - Dk[1];
 
 4266                Ddk[2] = dk[2] - Dk[2];
 
 4271    if (!discretised_field.is_null()) {
 
 4272        discretised_field.resize(number_of_dof);
 
 4273        discretised_field.set_zero();
 
 4276    if (!dof_numbering.is_null()) { get_dof_numbering(dof_numbering); }
 
 4279    for (
int l = 0; l != nb_gauss; ++l) {
 
 4280        const OnePoint& N = gauss_point[l];
 
 4284              b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin, e_end);
 
 4299        for (
int k = 0; k != nb_nodes; ++k) {
 
 4300            for (
int j = 0; j != 3; ++j) {
 
 4301                Ar[j] += N.dN[0][k] * R[k][j];
 
 4302                As[j] += N.dN[1][k] * R[k][j];
 
 4303                At[j] += N.N[k] * D[k][j];
 
 4304                Atr[j] += N.dN[0][k] * D[k][j];
 
 4305                Ats[j] += N.dN[1][k] * D[k][j];
 
 4306                ur[j] += N.dN[0][k] * u_dof[k][j];
 
 4307                us[j] += N.dN[1][k] * u_dof[k][j];
 
 4308                ut[j] += N.N[k] * Dd[k][j];
 
 4309                utr[j] += N.dN[0][k] * Dd[k][j];
 
 4310                uts[j] += N.dN[1][k] * Dd[k][j];
 
 4314        const double omega = (face == 4 ? -0.5 : (face == 5 ? +0.5 : 0.0));
 
 4318        TP d_coor_d_lcoor[2][3];
 
 4319        TP d_coor_deformed_d_lcoor[2][3];
 
 4320        for (
int j = 0; j != 3; ++j) {
 
 4321            d_coor_d_lcoor[0][j] = Ar[j] + omega * Atr[j];
 
 4322            d_coor_d_lcoor[1][j] = As[j] + omega * Ats[j];
 
 4323            d_coor_deformed_d_lcoor[0][j] = d_coor_d_lcoor[0][j] + ur[j] + omega * utr[j];
 
 4324            d_coor_deformed_d_lcoor[1][j] = d_coor_d_lcoor[1][j] + us[j] + omega * uts[j];
 
 4330            area = 
norm_3(normal) * gauss_point[l].weight;
 
 4338        const double xi[3] = {N.IntCoor[0], N.IntCoor[1], omega / (e_end - e_begin) * 2};
 
 4342        for (
int i = 0; i != nb_nodes; ++i) {
 
 4343            x_0[0] += N.N[i] * (R[i][0] + omega * D[i][0]);
 
 4344            x_0[1] += N.N[i] * (R[i][1] + omega * D[i][1]);
 
 4345            x_0[2] += N.N[i] * (R[i][2] + omega * D[i][2]);
 
 4349        TP x_d[3] = {x_0[0], x_0[1], x_0[2]};
 
 4350        for (
int i = 0; i != nb_nodes; ++i) {
 
 4351            x_d[0] += N.N[i] * (u_dof[i][0] + omega * (d[i][0] - D[i][0]));
 
 4352            x_d[1] += N.N[i] * (u_dof[i][1] + omega * (d[i][1] - D[i][1]));
 
 4353            x_d[2] += N.N[i] * (u_dof[i][2] + omega * (d[i][2] - D[i][2]));
 
 4357        b2linalg::Vector<TP, b2linalg::Vdense> value;
 
 4359              b2linalg::Vector<double, b2linalg::Vdense_constref>(3, xi),
 
 4360              b2linalg::Vector<TP, b2linalg::Vdense_constref>(3, x_0),
 
 4361              b2linalg::Vector<TP, b2linalg::Vdense_constref>(3, x_d),
 
 4362              b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>(3, 2, &d_coor_d_lcoor[0][0]),
 
 4363              b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>(
 
 4364                    3, 2, &d_coor_deformed_d_lcoor[0][0]),
 
 4375                const TP radius = 1.;
 
 4376                const TP sina = -y / radius;
 
 4377                const TP cosa = 1 - z / radius;
 
 4378                const TP cos2a = cosa * cosa - sina * sina;
 
 4379                const TP scale = cos2a;
 
 4382                const TP pi = std::acos(TP(-1));
 
 4384                    const TP ca = 0.8660254037844387;
 
 4389                const TP scale = std::sin(pi * x) * std::sin(pi * y);
 
 4394        TP Nl[max_nb_dof][3];
 
 4395        std::fill_n(Nl[0], number_of_dof * 3, 0);
 
 4397        for (
int k = 0; k != nb_nodes; ++k) {
 
 4398            const TP Nlk = N.N[k];
 
 4399            Nl[dofi][0] = Nl[dofi + 1][1] = Nl[dofi + 2][2] = Nlk;
 
 4401            const TP Nlk_omega = Nlk * omega;
 
 4402            if (nodes_type_6_dof[k]) {
 
 4404                Nl[dofi][1] = Nlk_omega * -D[k][2];
 
 4405                Nl[dofi][2] = Nlk_omega * D[k][1];
 
 4407                Nl[++dofi][0] = Nlk_omega * D[k][2];
 
 4409                Nl[dofi][2] = Nlk_omega * -D[k][0];
 
 4411                Nl[++dofi][0] = Nlk_omega * -D[k][1];
 
 4412                Nl[dofi][1] = Nlk_omega * D[k][0];
 
 4416                    const TP* 
const Dk = Dr[k][0];
 
 4417                    Nl[dofi][0] = Nlk_omega * (D[k][2] * Dk[1] - D[k][1] * Dk[2]);
 
 4418                    Nl[dofi][1] = Nlk_omega * (-D[k][2] * Dk[0] + D[k][0] * Dk[2]);
 
 4419                    Nl[dofi][2] = Nlk_omega * (D[k][1] * Dk[0] - D[k][0] * Dk[1]);
 
 4422                    const TP* 
const Dk = Dr[k][1];
 
 4423                    Nl[++dofi][0] = Nlk_omega * (D[k][2] * Dk[1] - D[k][1] * Dk[2]);
 
 4424                    Nl[dofi][1] = Nlk_omega * (-D[k][2] * Dk[0] + D[k][0] * Dk[2]);
 
 4425                    Nl[dofi][2] = Nlk_omega * (D[k][1] * Dk[0] - D[k][0] * Dk[1]);
 
 4432              'T', 3, number_of_dof, area, Nl[0], 3, &value[0], 1, 1.0, &discretised_field[0], 1);
 
 4437      typename TP, 
typename SHAPE, 
typename SHAPE_IP, 
bool SHAPE_IP_OOP_SAME,
 
 4438      typename MITC_TYING_POINT, 
typename DIRECTOR_ESTIMATION, 
bool incompatible_mode, 
bool mitc_gs,
 
 4439      int nb_gauss, 
typename INTEGRATION, 
char NAME_SUFFIX>
 
 4440void b2000::ShellMITC<
 
 4441      TP, SHAPE, SHAPE_IP, SHAPE_IP_OOP_SAME, MITC_TYING_POINT, DIRECTOR_ESTIMATION,
 
 4442      incompatible_mode, mitc_gs, nb_gauss, INTEGRATION, NAME_SUFFIX>::
 
 4443      field_volume_integration(
 
 4444            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dof,
 
 4445            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdot,
 
 4446            const b2linalg::Vector<TP, b2linalg::Vdense_constref>& dofdotdot,
 
 4447            const Field<TP>& field, b2linalg::Index& dof_numbering,
 
 4448            b2linalg::Vector<TP, b2linalg::Vdense>& discretised_field,
 
 4449            b2linalg::Matrix<TP, b2linalg::Mpacked>& d_discretised_field_d_dof) {
 
 4458#ifdef TT_INTEGRATION_AT_GAUSS_POINT 
 4459#if TT_INTEGRATION_AT_GAUSS_POINT == 2 
 4460    const double sqrt_inv_3 = std::sqrt(1 / 3.);
 
 4461    static const int nb_tt_integration = 2;
 
 4462    static const double n_omega[2] = {-sqrt_inv_3, sqrt_inv_3};
 
 4463    static const double n_weight[2] = {0.5, 0.5};
 
 4464#elif TT_INTEGRATION_AT_GAUSS_POINT == 3 
 4465    static const int nb_tt_integration = 3;
 
 4466    static const double n_omega[3] = {-1, 0, 1};
 
 4467    static const double n_weight[3] = {1. / 6, 4. / 6, 1. / 6};
 
 4470    static const int nb_tt_integration = 2;
 
 4471    static const double n_omega[2] = {-1, 1};
 
 4472    static const double n_weight[2] = {0.5, 0.5};
 
 4476    TP D[nb_nodes][3] = {};
 
 4477    TP Dd[nb_nodes][3] = {};
 
 4480    TP Dr[nb_nodes][2][3] = {};
 
 4483    TP R[nb_nodes][6] = {};
 
 4485    int number_of_dof = incompatible_mode ? 4 : 0;
 
 4487    for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 4488        if (nodes_type_6_dof[k]) {
 
 4490            dof6_node_type::get_coor_s(R[k], nodes[k]);
 
 4493            dof5_node_type::get_coor_s(R[k], nodes[k]);
 
 4497    for (
int k = nb_nodes_ip; k != nb_nodes; ++k) {
 
 4499        dof2_node_type::get_coor_s(R[k], nodes[k]);
 
 4503    DIRECTOR_ESTIMATION::compute_normal(R, D);
 
 4504    for (
int k = 0; k != nb_nodes_ip; ++k) {
 
 4505        if (!nodes_type_6_dof[k]
 
 4506            && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
 4507            if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
 4510            std::copy(R[k] + 3, R[k] + 6, D[k]);
 
 4515    for (
int k = 0; k != nb_nodes; ++k) {
 
 4516        if (!nodes_type_6_dof[k]) {
 
 4517            const TP e2[3] = {0, 1, 0};
 
 4520                const TP e2[3] = {1, 0, 0};
 
 4530    TP u_dof[nb_nodes][3];
 
 4533    TP w_dof[nb_nodes][3];
 
 4539    if (dof.is_null()) {
 
 4540        std::fill_n(u_dof[0], nb_nodes * 3, 0);
 
 4541        std::fill_n(w_dof[0], nb_nodes * 3, 0);
 
 4542        if (incompatible_mode) { std::fill_n(inc_dof, 4, 0); }
 
 4544        const TP* alpha = &dof[0];
 
 4545        for (
int k = 0; k != nb_nodes; ++k) {
 
 4546            if (nodes_type_6_dof[k]) {
 
 4547                std::copy(alpha, alpha + 3, u_dof[k]);
 
 4549                std::copy(alpha, alpha + 3, w_dof[k]);
 
 4552                std::copy(alpha, alpha + 3, u_dof[k]);
 
 4553                w_dof[k][0] = Dr[k][0][0] * alpha[3] + Dr[k][1][0] * alpha[4];
 
 4554                w_dof[k][1] = Dr[k][0][1] * alpha[3] + Dr[k][1][1] * alpha[4];
 
 4555                w_dof[k][2] = Dr[k][0][2] * alpha[3] + Dr[k][1][2] * alpha[4];
 
 4566    TP norm_w[nb_nodes];
 
 4570    for (
int k = 0; k != nb_nodes; ++k) {
 
 4571        const TP* w = w_dof[k];
 
 4572        const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
 4573        norm_w[k] = std::sqrt(ww);
 
 4576            const TP w0_2 = w[0] * w[0];
 
 4577            const TP w1_2 = w[1] * w[1];
 
 4578            const TP w2_2 = w[2] * w[2];
 
 4579            O2[0] = -w2_2 - w1_2;
 
 4580            O2[1] = w[0] * w[1];
 
 4581            O2[2] = w[0] * w[2];
 
 4582            O2[3] = -w2_2 - w0_2;
 
 4583            O2[4] = w[1] * w[2];
 
 4584            O2[5] = -w1_2 - w0_2;
 
 4593            s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
 4594            c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
 4596            s = std::sin(norm_w[k]) / norm_w[k];
 
 4597            c = std::sin(norm_w[k] * 0.5) / norm_w[k];
 
 4601            TP* 
const dk = d[k];
 
 4602            const TP* 
const Dk = D[k];
 
 4604                TP* 
const Ddk = Dd[k];
 
 4605                Ddk[0] = Dk[2] * w_dof[k][1] - Dk[1] * w_dof[k][2];
 
 4606                Ddk[1] = -Dk[2] * w_dof[k][0] + Dk[0] * w_dof[k][2];
 
 4607                Ddk[2] = Dk[1] * w_dof[k][0] - Dk[0] * w_dof[k][1];
 
 4608                dk[0] = Ddk[0] + Dk[0];
 
 4609                dk[1] = Ddk[1] + Dk[1];
 
 4610                dk[2] = Ddk[2] + Dk[2];
 
 4612                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 4613                        + (s * w[1] + c * O2[2]) * Dk[2];
 
 4614                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 4615                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
 4616                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 4617                        + (1 + c * O2[5]) * Dk[2];
 
 4618                TP* 
const Ddk = Dd[k];
 
 4619                Ddk[0] = dk[0] - Dk[0];
 
 4620                Ddk[1] = dk[1] - Dk[1];
 
 4621                Ddk[2] = dk[2] - Dk[2];
 
 4626    if (!discretised_field.is_null()) {
 
 4627        discretised_field.resize(number_of_dof);
 
 4628        discretised_field.set_zero();
 
 4631    if (!dof_numbering.is_null()) { get_dof_numbering(dof_numbering); }
 
 4634    for (
int l = 0; l != nb_gauss; ++l) {
 
 4635        const OnePoint& N = gauss_point[l];
 
 4638        const int number_of_layer = 
property->number_of_layer();
 
 4640              b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin, e_end);
 
 4655        for (
int k = 0; k != nb_nodes; ++k) {
 
 4656            for (
int j = 0; j != 3; ++j) {
 
 4657                Ar[j] += N.dN[0][k] * R[k][j];
 
 4658                As[j] += N.dN[1][k] * R[k][j];
 
 4659                At[j] += N.N[k] * D[k][j];
 
 4660                Atr[j] += N.dN[0][k] * D[k][j];
 
 4661                Ats[j] += N.dN[1][k] * D[k][j];
 
 4662                ur[j] += N.dN[0][k] * u_dof[k][j];
 
 4663                us[j] += N.dN[1][k] * u_dof[k][j];
 
 4664                ut[j] += N.N[k] * Dd[k][j];
 
 4665                utr[j] += N.dN[0][k] * Dd[k][j];
 
 4666                uts[j] += N.dN[1][k] * Dd[k][j];
 
 4670        for (
int layer = 0; layer != number_of_layer; ++layer) {
 
 4674                  b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
 
 4676            for (
int lt = 0; lt != nb_tt_integration; ++lt) {
 
 4677                const TP omega = l_begin + 0.5 * (1 + n_omega[lt]) * (l_end - l_begin);
 
 4680                for (
int j = 0; j != 3; ++j) {
 
 4681                    G_0[0][j] = Ar[j] + omega * Atr[j];
 
 4682                    G_0[1][j] = As[j] + omega * Ats[j];
 
 4684                    G_d[0][j] = ur[j] + omega * utr[j];
 
 4685                    G_d[1][j] = us[j] + omega * uts[j];
 
 4690                const double xi[3] = {N.IntCoor[0], N.IntCoor[1], omega / (e_end - e_begin) * 2};
 
 4694                for (
int i = 0; i != nb_nodes; ++i) {
 
 4695                    x_0[0] += N.N[i] * (R[i][0] + omega * D[i][0]);
 
 4696                    x_0[1] += N.N[i] * (R[i][1] + omega * D[i][1]);
 
 4697                    x_0[2] += N.N[i] * (R[i][2] + omega * D[i][2]);
 
 4701                TP x_d[3] = {x_0[0], x_0[1], x_0[2]};
 
 4702                for (
int i = 0; i != nb_nodes; ++i) {
 
 4703                    x_d[0] += N.N[i] * (u_dof[i][0] + omega * (d[i][0] - D[i][0]));
 
 4704                    x_d[1] += N.N[i] * (u_dof[i][1] + omega * (d[i][1] - D[i][1]));
 
 4705                    x_d[2] += N.N[i] * (u_dof[i][2] + omega * (d[i][2] - D[i][2]));
 
 4709                b2linalg::Vector<TP, b2linalg::Vdense> value;
 
 4711                      b2linalg::Vector<double, b2linalg::Vdense_constref>(3, xi),
 
 4712                      b2linalg::Vector<TP, b2linalg::Vdense_constref>(3, x_0),
 
 4713                      b2linalg::Vector<TP, b2linalg::Vdense_constref>(3, x_d),
 
 4714                      b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>(3, 3, &G_0[0][0]),
 
 4715                      b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>(3, 3, &G_d[0][0]), value);
 
 4719                if (field.get_multiply_with_density()) {
 
 4720                    mult = 
property->get_density(layer);
 
 4721                    add = (layer == 0 && lt == 1 ? non_structural_mass : TP(0));
 
 4726                const TP volume = n_weight[lt] * (l_end - l_begin) * det_G * gauss_point[l].weight;
 
 4729                TP Nl[max_nb_dof][3];
 
 4730                std::fill_n(Nl[0], number_of_dof * 3, 0);
 
 4732                for (
int k = 0; k != nb_nodes; ++k) {
 
 4733                    const TP Nlk = N.N[k];
 
 4734                    Nl[dofi][0] = Nl[dofi + 1][1] = Nl[dofi + 2][2] = Nlk;
 
 4736                    const TP Nlk_omega = Nlk * omega;
 
 4737                    if (nodes_type_6_dof[k]) {
 
 4739                        Nl[dofi][1] = Nlk_omega * -D[k][2];
 
 4740                        Nl[dofi][2] = Nlk_omega * D[k][1];
 
 4742                        Nl[++dofi][0] = Nlk_omega * D[k][2];
 
 4744                        Nl[dofi][2] = Nlk_omega * -D[k][0];
 
 4746                        Nl[++dofi][0] = Nlk_omega * -D[k][1];
 
 4747                        Nl[dofi][1] = Nlk_omega * D[k][0];
 
 4751                            const TP* 
const Dk = Dr[k][0];
 
 4752                            Nl[dofi][0] = Nlk_omega * (D[k][2] * Dk[1] - D[k][1] * Dk[2]);
 
 4753                            Nl[dofi][1] = Nlk_omega * (-D[k][2] * Dk[0] + D[k][0] * Dk[2]);
 
 4754                            Nl[dofi][2] = Nlk_omega * (D[k][1] * Dk[0] - D[k][0] * Dk[1]);
 
 4757                            const TP* 
const Dk = Dr[k][1];
 
 4758                            Nl[++dofi][0] = Nlk_omega * (D[k][2] * Dk[1] - D[k][1] * Dk[2]);
 
 4759                            Nl[dofi][1] = Nlk_omega * (-D[k][2] * Dk[0] + D[k][0] * Dk[2]);
 
 4760                            Nl[dofi][2] = Nlk_omega * (D[k][1] * Dk[0] - D[k][0] * Dk[1]);
 
 4767                      'T', 3, number_of_dof, (volume * mult + area * add), Nl[0], 3, &value[0], 1,
 
 4768                      1.0, &discretised_field[0], 1);
 
 4778      typename TP, 
typename SHAPE, 
typename DIRECTOR_ESTIMATION, 
typename EDGE, 
bool FREEZ = 
false>
 
 4781    using dof5_node_type = node::HNode<
 
 4783                coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
 
 4784                coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
 
 4786                coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
 
 4787                coordof::DTrans<coordof::RX, coordof::RY>>>;
 
 4789    using dof6_node_type = node::HNode<
 
 4791                coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
 
 4792                coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
 
 4794                coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
 
 4795                coordof::DTrans<coordof::RX, coordof::RY, coordof::RZ>>>;
 
 4797    using volume_node_type = node::HNode<
 
 4798          coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
 
 4799          coordof::Dof<coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>>>;
 
 4801    using type_t = ObjectTypeComplete<ShellMITC_Volume_Coupling, typename TypedElement<TP>::type_t>;
 
 4803    void set_nodes(std::pair<int, Node* const*> nodes_) {
 
 4804        if (nodes_.first != nb_nodes + 1) {
 
 4805            Exception() << 
"This element has " << nb_nodes + 1 << 
" nodes." << 
THROW;
 
 4807        for (
int i = 0; i != nb_nodes; ++i) {
 
 4808            nodes[i] = nodes_.second[i];
 
 4809            if (!dof5_node_type::is_dof_subset_s(nodes[i])) {
 
 4812                            << 
" because it does not have the dofs UX,UY,UZ,RX,RY." << 
THROW;
 
 4814            nodes_type_6_dof[i] = dof6_node_type::is_dof_subset_s(nodes[i]);
 
 4816        nodes[nb_nodes] = nodes_.second[nb_nodes];
 
 4817        if (!volume_node_type::is_dof_subset_s(nodes[nb_nodes])) {
 
 4819                        << 
" cannot accept as last node " << nodes[nb_nodes]->get_object_name()
 
 4820                        << 
" because it does " 
 4821                           "not have the dofs UX,UY,UZ." 
 4826    std::pair<int, Node* const*> get_nodes()
 const {
 
 4827        return std::pair<int, Node* const*>(nb_nodes + 1, nodes);
 
 4830    void get_dof_numbering(b2linalg::Index& dof_numbering) {
 
 4832        for (
int k = 0; k != EDGE::nb_nodes_of; ++k) {
 
 4833            s += nodes_type_6_dof[EDGE::nodes_of_id[k]] ? 6 : 5;
 
 4835        dof_numbering.resize(s + 3);
 
 4836        size_t* ptr = &dof_numbering[0];
 
 4837        for (
int k = 0; k != EDGE::nb_nodes_of; ++k) {
 
 4838            if (nodes_type_6_dof[EDGE::nodes_of_id[k]]) {
 
 4839                ptr = dof6_node_type::get_global_dof_numbering_s(ptr, nodes[EDGE::nodes_of_id[k]]);
 
 4841                ptr = dof5_node_type::get_global_dof_numbering_s(ptr, nodes[EDGE::nodes_of_id[k]]);
 
 4844        ptr = volume_node_type::get_global_dof_numbering_s(ptr, nodes[nb_nodes]);
 
 4847    void set_property(ElementProperty* property_) { 
property = property_; }
 
 4849    const ElementProperty* get_property()
 const { 
return property; }
 
 4851    void init(Model& model) {}
 
 4881    std::pair<int, Element::VariableInfo> get_constraint_info() {
 
 4882        return std::pair<int, Element::VariableInfo>(FREEZ ? 2 : 3, Element::
nonlinear);
 
 4885    void get_constraint(
 
 4886          Model& model, 
const bool linear, 
double time,
 
 4887          const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof,
 
 4888          b2linalg::Index& dof_numbering, b2linalg::Vector<TP, b2linalg::Vdense>& constraint,
 
 4889          b2linalg::Matrix<TP, b2linalg::Mcompressed_col>& trans_d_constraint_d_dof,
 
 4890          b2linalg::Vector<TP, b2linalg::Vdense>& d_constraint_d_time);
 
 4895    static const int nb_nodes = SHAPE::num_nodes;
 
 4896    ElementProperty* property;
 
 4897    Node* nodes[nb_nodes + 1];
 
 4901    std::bitset<nb_nodes> nodes_type_6_dof;
 
 4909          const TP a_coor[3], 
const TP b_coor[3], 
const TP c_coor[3], 
const TP d_coor[3], TP& va,
 
 4913        for (i = 0; i != 3; ++i) {
 
 4914            coef[0] += 4 * c_coor[i] * c_coor[i] - 4 * b_coor[i] * c_coor[i]
 
 4915                       - 4 * a_coor[i] * c_coor[i] + b_coor[i] * b_coor[i]
 
 4916                       + 2 * a_coor[i] * b_coor[i] + a_coor[i] * a_coor[i];
 
 4917            coef[1] += -3 * b_coor[i] * c_coor[i] + 3 * a_coor[i] * c_coor[i]
 
 4918                       + 1.5 * b_coor[i] * b_coor[i] - 1.5 * a_coor[i] * a_coor[i];
 
 4919            coef[2] += 4 * c_coor[i] * d_coor[i] - 2 * b_coor[i] * d_coor[i]
 
 4920                       - 2 * a_coor[i] * d_coor[i] - 4 * c_coor[i] * c_coor[i]
 
 4921                       + 2 * b_coor[i] * c_coor[i] + 2 * a_coor[i] * c_coor[i]
 
 4922                       + 0.5 * b_coor[i] * b_coor[i] - a_coor[i] * b_coor[i]
 
 4923                       + 0.5 * a_coor[i] * a_coor[i];
 
 4924            coef[3] += -b_coor[i] * d_coor[i] + a_coor[i] * d_coor[i] + b_coor[i] * c_coor[i]
 
 4925                       - a_coor[i] * c_coor[i];
 
 4934            if (r < -1) { r = -1; }
 
 4935            if (r > 1) { r = 1; }
 
 4938            TP l_border = 1e100;
 
 4939            for (j = 0; j != nsol; ++j) {
 
 4941                TP l_border_tmp = 0;
 
 4946                    l_border_tmp = -c - 1;
 
 4948                    l_border_tmp = c - 1;
 
 4951                coefd[0] = 0.5 * (c2 - c);
 
 4952                coefd[1] = 0.5 * (c2 + c);
 
 4954                for (i = 0; i != 3; ++i) {
 
 4955                    TP lr_tmp = a_coor[i] * coefd[0] + b_coor[i] * coefd[1] + c_coor[i] * coefd[2];
 
 4956                    lr2_tmp += lr_tmp * lr_tmp;
 
 4958                if (l_border_tmp < l_border && lr2_tmp < lr2) {
 
 4960                    l_border = l_border_tmp;
 
 4966        va = 0.5 * (r2 - r);
 
 4967        vb = 0.5 * (r2 + r);
 
 4975          const TP a_coor[3], 
const TP b_coor[3], 
const TP c_coor[3], 
const TP pos, TP tangent[3]) {
 
 4976        for (
int i = 0; i != 3; ++i) {
 
 4978                  (a_coor[i] + b_coor[i] - 2 * c_coor[i]) * pos + (b_coor[i] - a_coor[i]) / 2;
 
 4984template <
typename TP, 
typename SHAPE, 
typename DIRECTOR_ESTIMATION, 
typename EDGE, 
bool FREEZ>
 
 4985typename ShellMITC_Volume_Coupling<TP, SHAPE, DIRECTOR_ESTIMATION, EDGE, FREEZ>::type_t
 
 4986      ShellMITC_Volume_Coupling<TP, SHAPE, DIRECTOR_ESTIMATION, EDGE, FREEZ>::type(
 
 4987            "SSC." + std::string(SHAPE::name) + 
".S.MITC" + (FREEZ ? 
".TF" : 
"")
 
 4988                  + (std::is_same<TP, 
b2000::csda<double>>::value ? 
".CSDA" : 
""),
 
 4989            "", StringList(), element::module, &TypedElement<TP>::type);
 
 4992template <
typename TP, 
typename SHAPE, 
typename DIRECTOR_ESTIMATION, 
typename EDGE, 
bool FREEZ>
 
 4993void b2000::ShellMITC_Volume_Coupling<TP, SHAPE, DIRECTOR_ESTIMATION, EDGE, FREEZ>::get_constraint(
 
 4994      Model& model, 
const bool linear, 
double time,
 
 4995      const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof,
 
 4996      b2linalg::Index& dof_numbering, b2linalg::Vector<TP, b2linalg::Vdense>& constraint,
 
 4997      b2linalg::Matrix<TP, b2linalg::Mcompressed_col>& trans_d_constraint_d_dof,
 
 4998      b2linalg::Vector<TP, b2linalg::Vdense>& d_constraint_d_time) {
 
 4999    get_dof_numbering(dof_numbering);
 
 5001    const int nb_nodes_of = EDGE::nb_nodes_of;
 
 5002    const int* nodes_of_id = EDGE::nodes_of_id;
 
 5005    TP D[nb_nodes][3] = {};
 
 5008    TP Dr[nb_nodes_of][2][3] = {};
 
 5011    TP R[nb_nodes][6] = {};
 
 5013    for (
int k = 0; k != nb_nodes; ++k) {
 
 5014        if (nodes_type_6_dof[k]) {
 
 5015            dof6_node_type::get_coor_s(R[k], nodes[k]);
 
 5017            dof5_node_type::get_coor_s(R[k], nodes[k]);
 
 5021    int number_of_dof = 0;
 
 5022    for (
int k = 0; k != nb_nodes_of; ++k) {
 
 5023        if (nodes_type_6_dof[nodes_of_id[k]]) {
 
 5029    DIRECTOR_ESTIMATION::compute_normal(R, D);
 
 5030    for (
int k = 0; k != nb_nodes; ++k) {
 
 5031        if (!nodes_type_6_dof[k]
 
 5032            && (R[k][3] * R[k][3] + R[k][4] * R[k][4] + R[k][5] * R[k][5] > 0.5)) {
 
 5033            if (
dot_3(R[k] + 3, D[k]) < 0) {  
 
 5036            std::copy(R[k] + 3, R[k] + 6, D[k]);
 
 5041    for (
int k = 0; k != nb_nodes_of; ++k) {
 
 5042        const int kk = nodes_of_id[k];
 
 5044            std::copy(D[kk], D[kk + 1], D[k]);
 
 5045            std::copy(R[kk], R[kk + 1], R[k]);
 
 5049    for (
int kk = 0; kk != nb_nodes_of; ++kk) {
 
 5050        if (!nodes_type_6_dof[nodes_of_id[kk]]) {
 
 5051            TP e2[3] = {0, 1, 0};
 
 5054                TP e2[3] = {1, 0, 0};
 
 5064    TP u_dof[nb_nodes_of][3];
 
 5067    TP w_dof[nb_nodes_of][3];
 
 5070    if (dof.size2() == 0) {
 
 5071        std::fill_n(u_dof[0], nb_nodes_of * 3, 0);
 
 5072        std::fill_n(w_dof[0], nb_nodes_of * 3, 0);
 
 5075        for (
int k = 0; k != nb_nodes_of; ++k) {
 
 5076            const TP* alpha = &dof[0][dof_numbering[o]];
 
 5077            if (nodes_type_6_dof[k]) {
 
 5078                std::copy(alpha, alpha + 3, u_dof[k]);
 
 5079                std::copy(alpha + 3, alpha + 6, w_dof[k]);
 
 5082                std::copy(alpha, alpha + 3, u_dof[k]);
 
 5083                w_dof[k][0] = Dr[k][0][0] * alpha[3] + Dr[k][1][0] * alpha[4];
 
 5084                w_dof[k][1] = Dr[k][0][1] * alpha[3] + Dr[k][1][1] * alpha[4];
 
 5085                w_dof[k][2] = Dr[k][0][2] * alpha[3] + Dr[k][1][2] * alpha[4];
 
 5093        TP Dn[2][nb_nodes_of][3];
 
 5094        if (nb_nodes_of == 2) {
 
 5095            sub_3(R[0], R[1], Dn[0][0]);
 
 5097            std::copy(Dn[0][0], Dn[0][1], Dn[0][1]);
 
 5098        } 
else if (nb_nodes_of == 3) {
 
 5099            get_tangent(R[0], R[1], R[2], -1, Dn[0][0]);
 
 5100            get_tangent(R[0], R[1], R[2], 1, Dn[0][1]);
 
 5101            get_tangent(R[0], R[1], R[2], 0, Dn[0][2]);
 
 5107        TP d[2][nb_nodes_of][3];
 
 5111        TP T[2][nb_nodes_of][3][3];
 
 5115        for (
int kk = 0; kk != nb_nodes_of; ++kk) {
 
 5120            const TP* w = w_dof[kk];
 
 5121            const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
 5122            norm_w = std::sqrt(ww);
 
 5125                const TP w0_2 = w[0] * w[0];
 
 5126                const TP w1_2 = w[1] * w[1];
 
 5127                const TP w2_2 = w[2] * w[2];
 
 5128                O2[0] = -w2_2 - w1_2;
 
 5129                O2[1] = w[0] * w[1];
 
 5130                O2[2] = w[0] * w[2];
 
 5131                O2[3] = -w2_2 - w0_2;
 
 5132                O2[4] = w[1] * w[2];
 
 5133                O2[5] = -w1_2 - w0_2;
 
 5143                s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
 5144                c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
 5145                cc = 1.0 / 6 - ww / 120 + w4 / 5040 - w6 / 362880 + w8 / 39916800;
 
 5147                s = std::sin(norm_w) / norm_w;
 
 5148                c = std::sin(norm_w * 0.5) / norm_w;
 
 5150                cc = (norm_w - std::sin(norm_w)) / (ww * norm_w);
 
 5152            for (
int k = 0; k != 2; ++k) {
 
 5153                TP* 
const dk = d[k][kk];
 
 5154                const TP* 
const Dk = Dn[k][kk];
 
 5155                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 5156                        + (s * w[1] + c * O2[2]) * Dk[2];
 
 5157                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 5158                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
 5159                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 5160                        + (1 + c * O2[5]) * Dk[2];
 
 5163            if (trans_d_constraint_d_dof.is_null()) { 
continue; }
 
 5168            if (nodes_type_6_dof[nodes_of_id[kk]]) {
 
 5169                HT3k[0][0] = 1 + c * O2[0];
 
 5170                HT3k[1][0] = s * -w[2] + c * O2[1];
 
 5171                HT3k[2][0] = s * w[1] + c * O2[2];
 
 5172                HT3k[0][1] = s * w[2] + c * O2[1];
 
 5173                HT3k[1][1] = 1 + c * O2[3];
 
 5174                HT3k[2][1] = s * -w[0] + c * O2[4];
 
 5175                HT3k[0][2] = s * -w[1] + c * O2[2];
 
 5176                HT3k[1][2] = s * w[0] + c * O2[4];
 
 5177                HT3k[2][2] = 1 + c * O2[5];
 
 5180                    const TP* 
const Dk = Dr[kk][0];
 
 5181                    HT3k[0][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 5182                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
 5183                    HT3k[0][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 5184                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
 5185                    HT3k[0][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 5186                                 + (1 + c * O2[5]) * Dk[2];
 
 5189                    const TP* 
const Dk = Dr[kk][1];
 
 5190                    HT3k[1][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 5191                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
 5192                    HT3k[1][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 5193                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
 5194                    HT3k[1][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 5195                                 + (1 + c * O2[5]) * Dk[2];
 
 5198            for (
int k = 0; k != 2; ++k) {
 
 5199                const TP* 
const dk = d[k][kk];
 
 5200                T[k][kk][0][0] = dk[2] * HT3k[0][1] - dk[1] * HT3k[0][2];
 
 5201                T[k][kk][0][1] = -dk[2] * HT3k[0][0] + dk[0] * HT3k[0][2];
 
 5202                T[k][kk][0][2] = dk[1] * HT3k[0][0] - dk[0] * HT3k[0][1];
 
 5203                T[k][kk][1][0] = dk[2] * HT3k[1][1] - dk[1] * HT3k[1][2];
 
 5204                T[k][kk][1][1] = -dk[2] * HT3k[1][0] + dk[0] * HT3k[1][2];
 
 5205                T[k][kk][1][2] = dk[1] * HT3k[1][0] - dk[0] * HT3k[1][1];
 
 5206                if (nodes_type_6_dof[nodes_of_id[kk]]) {
 
 5207                    T[k][kk][2][0] = dk[2] * HT3k[2][1] - dk[1] * HT3k[2][2];
 
 5208                    T[k][kk][2][1] = -dk[2] * HT3k[2][0] + dk[0] * HT3k[2][2];
 
 5209                    T[k][kk][2][2] = dk[1] * HT3k[2][0] - dk[0] * HT3k[2][1];
 
 5215        volume_node_type::get_coor_s(c_coor, nodes[nb_nodes]);
 
 5218        if (nb_nodes_of == 2) {
 
 5221            for (
int i = 0; i != 3; ++i) {
 
 5222                const TP b_a = R[1][i] - R[0][i];
 
 5224                v[1] += b_a * (c_coor[i] - R[0][i]);
 
 5228        } 
else if (nb_nodes_of == 3) {
 
 5229            get_p_weight(R[0], R[1], R[2], c_coor, v[0], v[1], v[2]);
 
 5237            assert((
size_t)number_of_dof < dof_numbering.size());
 
 5238            const TP* c_disp = dof.is_null() ? 0 : &dof(dof_numbering[number_of_dof], 0);
 
 5239            for (
int i = 0; i != 3; ++i) {
 
 5240                cu[i] = (c_disp ? c_disp[i] : TP(0)) + c_coor[i];
 
 5241                cd[0][i] = cd[1][i] = 0;
 
 5242                for (
int k = 0; k != nb_nodes_of; ++k) {
 
 5243                    cu[i] -= (u_dof[k][i] + R[k][i]) * v[k];
 
 5244                    cd[0][i] += d[0][k][i] * v[k];
 
 5245                    cd[1][i] += d[1][k][i] * v[k];
 
 5250        if (!constraint.is_null()) {
 
 5251            constraint.resize(2);
 
 5252            if (dof.size2() == 0) {
 
 5253                constraint.set_zero();
 
 5255                for (
int k = 0; k != 2; ++k) {
 
 5257                    for (
int i = 0; i != 3; ++i) { constraint[k] += cu[i] * cd[k][i]; }
 
 5262        if (!trans_d_constraint_d_dof.is_null()) {
 
 5263            trans_d_constraint_d_dof.set_zero();
 
 5264            trans_d_constraint_d_dof.resize(number_of_dof + 3, 2, 2 * (number_of_dof + 3));
 
 5268            trans_d_constraint_d_dof.get_values(colind, rowind, value);
 
 5269            size_t* rowind_begin = rowind;
 
 5271            for (
int i = 0; i != 2; ++i) {
 
 5273                for (
int k = 0; k != nb_nodes_of; ++k) {
 
 5275                    for (
int j = 0; j != 3; ++j, ++pos) {
 
 5277                        *value++ = -v[k] * d[i][k][j];
 
 5278                        dd[j] = cu[j] * v[k];
 
 5280                    int s = nodes_type_6_dof[nodes_of_id[k]] ? 3 : 2;
 
 5281                    for (
int j = 0; j != s; ++j, ++pos) {
 
 5284                        for (
int jj = 0; jj != 3; ++jj) { tmp += dd[jj] * T[i][k][j][jj]; }
 
 5288                for (
int j = 0; j != 3; ++j, ++pos) {
 
 5290                    *value++ = cd[i][j];
 
 5292                *colind++ = rowind - rowind_begin;
 
 5298        TP d[nb_nodes_of][3];
 
 5302        TP T[nb_nodes_of][3][3];
 
 5306        for (
int kk = 0; kk != nb_nodes_of; ++kk) {
 
 5308            const TP* w = w_dof[kk];
 
 5309            const TP ww = w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
 
 5310            norm_w = std::sqrt(ww);
 
 5313                const TP w0_2 = w[0] * w[0];
 
 5314                const TP w1_2 = w[1] * w[1];
 
 5315                const TP w2_2 = w[2] * w[2];
 
 5316                O2[0] = -w2_2 - w1_2;
 
 5317                O2[1] = w[0] * w[1];
 
 5318                O2[2] = w[0] * w[2];
 
 5319                O2[3] = -w2_2 - w0_2;
 
 5320                O2[4] = w[1] * w[2];
 
 5321                O2[5] = -w1_2 - w0_2;
 
 5331                s = 1.0 - ww / 6 + w4 / 120 - w6 / 5040 + w8 / 362880 - w10 / 39916800;
 
 5332                c = 0.5 - ww / 24 + w4 / 720 - w6 / 40320 + w8 / 3628800 - w10 / 479001600;
 
 5333                cc = 1.0 / 6 - ww / 120 + w4 / 5040 - w6 / 362880 + w8 / 39916800;
 
 5335                s = std::sin(norm_w) / norm_w;
 
 5336                c = std::sin(norm_w * 0.5) / norm_w;
 
 5338                cc = (norm_w - std::sin(norm_w)) / (ww * norm_w);
 
 5341                TP* 
const dk = d[kk];
 
 5342                const TP* 
const Dk = D[kk];
 
 5343                dk[0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 5344                        + (s * w[1] + c * O2[2]) * Dk[2];
 
 5345                dk[1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 5346                        + (s * -w[0] + c * O2[4]) * Dk[2];
 
 5347                dk[2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 5348                        + (1 + c * O2[5]) * Dk[2];
 
 5351            if (trans_d_constraint_d_dof.is_null()) { 
continue; }
 
 5356            if (nodes_type_6_dof[nodes_of_id[kk]]) {
 
 5357                HT3k[0][0] = 1 + c * O2[0];
 
 5358                HT3k[1][0] = s * -w[2] + c * O2[1];
 
 5359                HT3k[2][0] = s * w[1] + c * O2[2];
 
 5360                HT3k[0][1] = s * w[2] + c * O2[1];
 
 5361                HT3k[1][1] = 1 + c * O2[3];
 
 5362                HT3k[2][1] = s * -w[0] + c * O2[4];
 
 5363                HT3k[0][2] = s * -w[1] + c * O2[2];
 
 5364                HT3k[1][2] = s * w[0] + c * O2[4];
 
 5365                HT3k[2][2] = 1 + c * O2[5];
 
 5368                    const TP* 
const Dk = Dr[kk][0];
 
 5369                    HT3k[0][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 5370                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
 5371                    HT3k[0][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 5372                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
 5373                    HT3k[0][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 5374                                 + (1 + c * O2[5]) * Dk[2];
 
 5377                    const TP* 
const Dk = Dr[kk][1];
 
 5378                    HT3k[1][0] = (1 + c * O2[0]) * Dk[0] + (s * -w[2] + c * O2[1]) * Dk[1]
 
 5379                                 + (s * w[1] + c * O2[2]) * Dk[2];
 
 5380                    HT3k[1][1] = (s * w[2] + c * O2[1]) * Dk[0] + (1 + c * O2[3]) * Dk[1]
 
 5381                                 + (s * -w[0] + c * O2[4]) * Dk[2];
 
 5382                    HT3k[1][2] = (s * -w[1] + c * O2[2]) * Dk[0] + (s * w[0] + c * O2[4]) * Dk[1]
 
 5383                                 + (1 + c * O2[5]) * Dk[2];
 
 5387                const TP* 
const dk = d[kk];
 
 5388                T[kk][0][0] = dk[2] * HT3k[0][1] - dk[1] * HT3k[0][2];
 
 5389                T[kk][0][1] = -dk[2] * HT3k[0][0] + dk[0] * HT3k[0][2];
 
 5390                T[kk][0][2] = dk[1] * HT3k[0][0] - dk[0] * HT3k[0][1];
 
 5391                T[kk][1][0] = dk[2] * HT3k[1][1] - dk[1] * HT3k[1][2];
 
 5392                T[kk][1][1] = -dk[2] * HT3k[1][0] + dk[0] * HT3k[1][2];
 
 5393                T[kk][1][2] = dk[1] * HT3k[1][0] - dk[0] * HT3k[1][1];
 
 5394                if (nodes_type_6_dof[nodes_of_id[kk]]) {
 
 5395                    T[kk][2][0] = dk[2] * HT3k[2][1] - dk[1] * HT3k[2][2];
 
 5396                    T[kk][2][1] = -dk[2] * HT3k[2][0] + dk[0] * HT3k[2][2];
 
 5397                    T[kk][2][2] = dk[1] * HT3k[2][0] - dk[0] * HT3k[2][1];
 
 5403        volume_node_type::get_coor_s(c_coor, nodes[nb_nodes]);
 
 5406        if (nb_nodes_of == 2) {
 
 5409            for (
int i = 0; i != 3; ++i) {
 
 5410                const TP b_a = R[1][i] - R[0][i];
 
 5412                v[1] += b_a * (c_coor[i] - R[0][i]);
 
 5416        } 
else if (nb_nodes_of == 3) {
 
 5417            get_p_weight(R[0], R[1], R[2], c_coor, v[0], v[1], v[2]);
 
 5423        for (
int i = 0; i != 3; ++i) {
 
 5426            for (
int j = 0; j != nb_nodes_of; ++j) {
 
 5427                p -= R[j][i] * v[j];
 
 5428                dir += D[j][i] * v[j];
 
 5433        if (!constraint.is_null()) {
 
 5434            constraint.resize(3);
 
 5435            if (dof.size2() == 0) {
 
 5436                constraint.set_zero();
 
 5438                const TP* c_disp = &dof(dof_numbering[number_of_dof], 0);
 
 5439                for (
int i = 0; i != 3; ++i) {
 
 5440                    TP ctmp = -c_disp[i];
 
 5441                    for (
int j = 0; j != nb_nodes_of; ++j) {
 
 5442                        ctmp += (u_dof[j][i] + (d[j][i] - D[j][i]) * vz) * v[j];
 
 5444                    constraint[i] = ctmp;
 
 5449        if (!trans_d_constraint_d_dof.is_null()) {
 
 5450            trans_d_constraint_d_dof.set_zero();
 
 5451            trans_d_constraint_d_dof.resize(
 
 5452                  number_of_dof + 3, 3, 3 * number_of_dof - 6 * nb_nodes_of + 3);
 
 5456            trans_d_constraint_d_dof.get_values(colind, rowind, value);
 
 5457            size_t* rowind_begin = rowind;
 
 5459            for (
int i = 0; i != 3; ++i) {
 
 5461                for (
int k = 0; k != nb_nodes_of; ++k) {
 
 5462                    *rowind++ = i + pos;
 
 5465                    int s = nodes_type_6_dof[nodes_of_id[k]] ? 3 : 2;
 
 5466                    for (
int j = 0; j != s; ++j) {
 
 5467                        *rowind++ = pos + j;
 
 5468                        *value++ = T[k][j][i] * vz * v[k];
 
 5472                *rowind++ = pos + i;
 
 5474                *colind++ = rowind - rowind_begin;
 
 5480#ifdef TT_INTEGRATION_AT_GAUSS_POINT 
 5481#undef TT_INTEGRATION_AT_GAUSS_POINT 
#define THROW
Definition b2exception.H:198
 
virtual void face_geom(const int face_id, const b2linalg::Vector< double, b2linalg::Vdense_constref > &internal_coor, b2linalg::Vector< double > &geom, b2linalg::Matrix< double, b2linalg::Mrectangle > &d_geom_d_icoor)
Definition b2element.H:876
 
virtual void edge_geom(const int edge_id, const double internal_coor, b2linalg::Vector< double > &geom, b2linalg::Vector< double > &d_geom_d_icoor)
Definition b2element.H:770
 
@ nonlinear
Definition b2element.H:619
 
@ linear
Definition b2element.H:615
 
const std::string & get_object_name() const override
Definition b2element.H:220
 
virtual void body_geom(const b2linalg::Vector< double, b2linalg::Vdense_constref > &internal_coor, b2linalg::Vector< double > &geom, b2linalg::Matrix< double, b2linalg::Mrectangle > &d_geom_d_icoor)
Definition b2element.H:963
 
Definition b2element.H:1054
 
virtual void get_value(Model &model, const bool linear, const EquilibriumSolution equilibrium_solution, const double time, const double delta_time, const b2linalg::Matrix< TP, b2linalg::Mrectangle_constref > &dof, GradientContainer *gradient_container, SolverHints *solver_hints, b2linalg::Index &dof_numbering, b2linalg::Vector< TP, b2linalg::Vdense > &value, const std::vector< bool > &d_value_d_dof_flags, std::vector< b2linalg::Matrix< TP, b2linalg::Mrectangle > > &d_value_d_dof, b2linalg::Vector< TP, b2linalg::Vdense > &d_value_d_time)
Compute the internal forces and their derivatives of the element.
Definition b2element.H:1619
 
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
 
T norm_outer_product_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:407
 
T norm_3(T a[3])
Definition b2tensor_calculus.H:364
 
void inner_product_3_3_NN(const T a[3][3], const T b[3][3], T c[3][3])
Definition b2tensor_calculus.H:808
 
T determinant_3_3(const T a[3][3])
Definition b2tensor_calculus.H:669
 
void sub_3(const T1 v1[3], const T2 v2[3], T3 dest[3])
Definition b2tensor_calculus.H:239
 
T normalise_3(T a[3])
Definition b2tensor_calculus.H:418
 
void outer_product_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:389
 
void solve_cubic_equation(const T a, const T b, const T c, const T d, int &nsol, T &x1, T &x2, T &x3)
Definition b2util.H:770
 
T invert_3_3(const T a[3][3], T b[3][3])
Definition b2tensor_calculus.H:699
 
T dot_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:328
 
void scale_3(T1 v[3], const T2 s)
Definition b2tensor_calculus.H:289
 
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314