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