b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_shell_mitc.H
1//------------------------------------------------------------------------
2// b2element_shell_mitc.H --
3//
4//
5// written by Mathias Doreille
6// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
7// Thomas Blome <thomas.blome@dlr.de>
8//
9// (c) 2007-2016,2017 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#ifndef B2ELEMENT_SHELL_MITC_H_
22#define B2ELEMENT_SHELL_MITC_H_
23
24#include <bitset>
25#include <cmath>
26#include <string>
27
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"
33#include "utils/b2exception.H"
34#include "utils/b2linear_algebra.H"
35#include "utils/b2rtable.H"
37
38// #define CHECK_TYING_INTERPOLATION 1
39#define TT_INTEGRATION_AT_GAUSS_POINT 3
40// #define TAYLOR_IM
41
42namespace b2000 {
43
44struct MITCTyingPointDef {
45 static const int nb_tying_point = 0;
46
47 struct TyingPoint {
48 double r;
49 double s;
50 bool comp_ip;
51 bool comp_oop;
52 };
53 static const TyingPoint tying_point[nb_tying_point];
54
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;
60
61 struct Interp {
62 int tp_pos;
63 int e_pos;
64 double weight;
65 };
66
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];
73 };
74
75 static void get_interpolation(double r, double s, TyingPointInterp& interp) {}
76};
77
78template <
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>
82class ShellMITC : public b2000::TypedElement<TP> {
83public:
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>>>;
87
88 using dof5_node_type = node::HNode<
89 coordof::Coor<
90 coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
91 coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
92 coordof::Dof<
93 coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
94 coordof::DTrans<coordof::RX, coordof::RY>>>;
95
96 using dof6_node_type = node::HNode<
97 coordof::Coor<
98 coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
99 coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
100 coordof::Dof<
101 coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
102 coordof::DTrans<coordof::RX, coordof::RY, coordof::RZ>>>;
103
104 using type_t = ObjectTypeComplete<ShellMITC, typename TypedElement<TP>::type_t>;
105
106 ShellMITC()
107 : incompatible_dof_index(0),
108 property(0),
109 non_structural_mass(0),
110 line_load_as_moment(false),
111 properties_list(0) {}
112
113 ~ShellMITC() {
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];
122 }
123 }
124 }
125 }
126 delete[] properties_list;
127 }
128 }
129
130 // The number of dof of the element (the incompatible dof)
131 int get_number_of_dof() const override {
132 if (incompatible_mode) {
133 return 4;
134 } else {
135 return 0;
136 }
137 }
138
139 size_t set_global_dof_numbering(size_t index) override {
140 incompatible_dof_index = index;
141 if (incompatible_mode) {
142 return index + 4;
143 } else {
144 return index;
145 }
146 }
147
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);
151 } else {
152 return std::pair<size_t, size_t>(0, 0);
153 }
154 }
155
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;
159 }
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])) {
163 Exception() << "MITC shell element " << this->get_object_name()
164 << " cannot accept node " << nodes[i]->get_object_name()
165 << " because it does not have the dofs UX,UY,UZ,RX,RY." << THROW;
166 }
167 nodes_type_6_dof[i] = dof6_node_type::is_dof_subset_s(nodes[i]);
168 }
169 for (int i = nb_nodes_ip; i != nb_nodes; ++i) {
170 if (!dof2_node_type::is_dof_subset_s(nodes[i])) {
171 Exception() << "MITC shell element " << this->get_object_name()
172 << " cannot accept node " << nodes[i]->get_object_name()
173 << " because it does not have the dofs RX,RY." << THROW;
174 }
175 nodes_type_6_dof[i] = false;
176 }
177 }
178
179 std::pair<int, Node* const*> get_nodes() const override {
180 return std::pair<int, Node* const*>(nb_nodes, nodes);
181 }
182
183 void set_property(ElementProperty* property_) override {
184 property = dynamic_cast<SolidMechanics25D<TP>*>(property_);
185 if (property == 0) {
186 Exception() << "Bad or missing element property type for shell "
187 "element "
188 << this->get_object_name() << " (forgot MID or PID element attribute?) "
189 << THROW;
190 }
191 }
192
193 const ElementProperty* get_property() const override { return property; }
194
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; }
200 }
201 }
202
203 void add_dof_field(Element::ListDofField& ldfield) override {
204 if (incompatible_mode) {
205 ldfield.add(
206 "DISPLACEMENT",
207 std::pair<size_t, size_t>(incompatible_dof_index, incompatible_dof_index + 4));
208 }
209 }
210
211 void init(Model& model) override {
212 if (property->linear() == SolidMechanics::yes) {
213 delete[] C_linear;
214 C_linear = 0;
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];
223 }
224 }
225 }
226 delete[] properties_list;
227 }
228 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));
235 }
236 }
237 }
238 } else {
239 properties_list = 0;
240 }
241
242 //
243 // Check for inverted element in the initial
244 // configuration, using double precision.
245 //
246 {
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]);
251 } else {
252 dof5_node_type::get_coor_s(R[k], nodes[k]);
253 }
254 }
255
256 const double eps = 1e-12;
257 bool inverted = false;
258
259 // Calculate surface normal at element center (not normalised).
260 double n_center[3];
261 {
262 double xi[3] = {};
263 b2linalg::Matrix<double> dN(2, nb_nodes);
264 SHAPE::eval_h_derived(xi, dN);
265
266 double Ar[3] = {};
267 double As[3] = {};
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];
272 }
273 }
274 outer_product_3(Ar, As, n_center);
275 if (normalise_3(n_center) < eps) { inverted = true; }
276 }
277
278 INTEGRATION integration;
279 integration.set_num(nb_gauss);
280 for (int l = 0; !inverted && l != nb_gauss; ++l) {
281 IP_ID ip_id;
282 integration.get_point(l, ip_id, gauss_point[l].weight);
283
284 b2linalg::Matrix<double> dN(2, nb_nodes);
285 SHAPE::eval_h_derived(ip_id.data(), dN);
286
287 // Calculate the surface normal at the integration
288 // point (not normalised).
289 double Ar[3] = {};
290 double As[3] = {};
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];
295 }
296 }
297 double n[3];
298 outer_product_3(Ar, As, n);
299 if (normalise_3(n) < eps) {
300 inverted = true;
301 break;
302 }
303
304 // Compare it with the surface normal at the center.
305 if (dot_3(n, n_center) < eps) {
306 inverted = true;
307 break;
308 }
309 }
310 if (inverted) {
311 Exception() << "The MITC shell element " << this->get_object_name()
312 << " is inverted." << THROW;
313 }
314 }
315 }
316
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]);
323 } else {
324 ptr = dof5_node_type::get_global_dof_numbering_s(ptr, nodes[k]);
325 }
326 }
327 for (int k = nb_nodes_ip; k != nb_nodes; ++k) {
328 ptr = dof2_node_type::get_global_dof_numbering_s(ptr, nodes[k]);
329 }
330 if (incompatible_mode) {
331 for (int i = 0; i != 4; ++i) { ptr[i] = incompatible_dof_index + i; }
332 }
333 }
334
335 const std::vector<Element::VariableInfo> get_value_info() const override {
336 return std::vector<Element::VariableInfo>(3, Element::nonlinear);
337 }
338
339 void get_value(
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;
348
349 int edge_field_order(const int edge_id, const std::string& field_name) override {
350 return SHAPE::edge_shape_type_0::order;
351 }
352
353 bool edge_field_linear_on_dof(const int edge_id, const std::string& field_name) override {
354 return true;
355 }
356
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 {
360 sub_nodes.clear();
361 sub_nodes.reserve(2);
362 sub_nodes.push_back(-1);
363 sub_nodes.push_back(1);
364 return 1;
365 }
366
367 void edge_geom(
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]]);
376 } else {
377 dof5_node_type::get_coor_s(node_coor[k], nodes[edge_node[k]]);
378 }
379 }
380 if (!geom.is_null()) {
381 b2linalg::Vector<double> shape(edge_shape::num_nodes);
382 edge_shape::eval_h(&internal_coor, shape);
383 geom.resize(3);
384 for (size_t i = 0; i != 3; ++i) {
385 TP tmp = 0;
386 for (int k = 0; k != edge_shape::num_nodes; ++k) {
387 tmp += shape[k] * node_coor[k][i];
388 }
389 geom[i] = tmp;
390 }
391 }
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) {
397 TP tmp = 0;
398 for (int k = 0; k != edge_shape::num_nodes; ++k) {
399 tmp += d_shape(0, k) * node_coor[k][i];
400 }
401 d_geom_d_icoor[i] = tmp;
402 }
403 }
404 }
405
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];
416 {
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;
421 } else {
422 ptr = dof5_node_type::get_global_dof_numbering_s(ptr, nodes[edge_node[k]]) - 2;
423 }
424 }
425 }
426 if (!value.is_null()) {
427 b2linalg::Vector<double> shape(edge_shape::num_nodes);
428 edge_shape::eval_h(&internal_coor, shape);
429 value.resize(3);
430 for (size_t i = 0; i != 3; ++i) {
431 TP tmp = 0;
432 for (int k = 0; k != edge_shape::num_nodes; ++k) {
433 tmp += shape[k] * dof(dn[k * 3 + i], 0);
434 }
435 value[i] = tmp;
436 }
437 }
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) {
443 TP tmp = 0;
444 for (int k = 0; k != edge_shape::num_nodes; ++k) {
445 tmp += d_shape(0, k) * dof(dn[k * 3 + i], 0);
446 }
447 d_value_d_icoor[i] = tmp;
448 }
449 }
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]; }
459 }
460 }
461 }
462
463 int face_field_order(const int face_id, const std::string& field_name) {
464 if (face_id > 4) {
465 return SHAPE::order;
466 } else {
467 return SHAPE::edge_shape_type_0::order + 1;
468 }
469 }
470
471 bool face_field_linear_on_dof(const int face_id, const std::string& field_name) override {
472 // due to the rotation dof, the displacement on the face is nonlinear
473 // if face_id is not 7 or 8
474 return face_id > 6;
475 }
476
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 {
481 // we split the rectangular face in two triangles
482 sub_nodes.resize(2, 4);
483 sub_nodes(0, 0) = -1;
484 sub_nodes(1, 0) = -1;
485 sub_nodes(0, 1) = 1;
486 sub_nodes(1, 1) = -1;
487 sub_nodes(0, 2) = 1;
488 sub_nodes(1, 2) = 1;
489 sub_nodes(0, 3) = -1;
490 sub_nodes(1, 3) = 1;
491 sub_faces.clear();
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);
496 }
497
498 void face_geom(
499 const int face_id,
500 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
501 b2linalg::Vector<TP>& geom, b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_geom_d_icoor) {
502 if (nb_nodes != nb_nodes_ip) { UnimplementedError() << THROW; }
503
504 typedef typename SHAPE::edge_shape_type_0 edge_shape;
505 TP R[nb_nodes][6] = {};
506 TP D[nb_nodes][3] = {};
507 TP e_begin, e_end;
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]);
511 } else {
512 dof5_node_type::get_coor_s(R[k], nodes[k]);
513 }
514 }
515 if (face_id < 7) {
516 // the surface is not the mid surface, we must also take into account the
517
518 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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) { // inversed normal
524 scale_3(R[k] + 3, -1.);
525 }
526 std::copy(R[k] + 3, R[k] + 6, D[k]);
527 }
528 }
529 b2linalg::Vector<double> N(nb_nodes);
530 SHAPE::eval_h(&internal_coor[0], N);
531 property->laminate(N, e_begin, e_end);
532 } else {
533 std::fill_n(D[0], 3 * nb_nodes, 0.0);
534 }
535
536 const int* map_node;
537 int nb_r_nodes;
538 TP z;
539 if (face_id > 4) {
540 if (face_id < 7) {
541 z = face_id == 5 ? e_begin : e_end;
542 map_node = SHAPE::face_node[face_id - 5];
543 } else {
544 z = 0;
545 map_node = SHAPE::face_node[face_id - 7];
546 }
547 nb_r_nodes = nb_nodes;
548 } else {
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);
552 }
553
554 if (!geom.is_null()) {
555 b2linalg::Vector<double> shape(nb_r_nodes);
556 if (face_id > 4) {
557 SHAPE::eval_h(&internal_coor[0], shape);
558 } else {
559 edge_shape::eval_h(&internal_coor[0], shape);
560 }
561 geom.resize(3);
562 for (size_t i = 0; i != 3; ++i) {
563 TP tmp = 0;
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]);
567 }
568 geom[i] = tmp;
569 }
570 }
571 if (!d_geom_d_icoor.is_null()) {
572 if (face_id > 4) {
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) {
578 TP tmp = 0;
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]);
582 }
583 d_geom_d_icoor(i, j) = tmp;
584 }
585 }
586 } else {
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) {
594 TP tmp0 = 0;
595 TP tmp1 = 0;
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];
600 }
601 d_geom_d_icoor(i, 0) = tmp0;
602 d_geom_d_icoor(i, 1) = tmp1 * z_scale;
603 }
604 }
605 }
606 }
607
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) {
616 if (nb_nodes != nb_nodes_ip) { UnimplementedError() << THROW; }
617
618 typedef typename SHAPE::edge_shape_type_0 edge_shape;
619 TP R[nb_nodes][6] = {};
620 TP D[nb_nodes][3] = {};
621
622 // The local referential for the nodes that have 5 dofs.
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]);
626 } else {
627 dof5_node_type::get_coor_s(R[k], nodes[k]);
628 }
629 }
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) { // inversed normal w.r.t. element
635 scale_3(R[k] + 3, -1.);
636 }
637 std::copy(R[k] + 3, R[k] + 6, D[k]);
638 }
639 }
640
641 TP z;
642 TP z_scale = 0;
643 {
644 TP e_begin, e_end;
645 if (face_id < 7) {
646 b2linalg::Vector<double> N(nb_nodes);
647 SHAPE::eval_h(&internal_coor[0], N);
648 property->laminate(N, e_begin, e_end);
649 if (face_id < 5) {
650 z = 0.5 * ((e_end - e_begin) * internal_coor[1] + e_end + e_begin);
651 z_scale = 0.5 * (e_end - e_begin);
652 } else {
653 z = face_id == 5 ? e_begin : e_end;
654 }
655 } else {
656 z = 0;
657 }
658 }
659
660 int nb_nodes_f;
661 const int* map_node;
662 if (face_id < 5) {
663 nb_nodes_f = edge_shape::num_nodes;
664 map_node = SHAPE::edge_node[face_id - 1];
665
666 } else {
667 nb_nodes_f = SHAPE::num_nodes;
668 if (face_id < 7) {
669 map_node = SHAPE::face_node[face_id - 5];
670 } else {
671 map_node = SHAPE::face_node[face_id - 7];
672 }
673 }
674 // update D and R to only contain the nodes of the face
675 if (face_id < 7) {
676 {
677 TP tmp[nb_nodes][6];
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]);
680 }
681 std::copy(tmp[0], tmp[edge_shape::num_nodes], R[0]);
682 }
683 {
684 TP tmp[nb_nodes][3];
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]);
687 }
688 std::copy(tmp[0], tmp[edge_shape::num_nodes], D[0]);
689 }
690 }
691
692 TP Dr[nb_nodes][2][3];
693
694 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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};
698 outer_product_3(e2, D[k], Dr[k][0]);
699 if (normalise_3(Dr[k][0]) < 0.01) {
700 const TP e2[3] = {1, 0, 0};
701 outer_product_3(e2, D[k], Dr[k][0]);
702 normalise_3(Dr[k][0]);
703 }
704 outer_product_3(D[k], Dr[k][0], Dr[k][1]);
705 normalise_3(Dr[k][1]);
706 }
707 }
708
709 // The 3 translations at each nodes
710 TP u_dof[nb_nodes][3];
711
712 // The 3 rotations at each nodes
713 TP w_dof[nb_nodes][3];
714
715 // The incompatible dof
716 TP inc_dof[4];
717
718 std::vector<size_t> dnvec(nb_nodes_f * 6 + 4);
719 size_t* dn{dnvec.data()};
720 size_t* ptr_dn = dn;
721
722 // convert the 5 dof to 6 dof
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; }
732 } else {
733 ptr_dn = dof5_node_type::get_global_dof_numbering_s(ptr_dn, nodes[kk]);
734 if (face_id > 6) { ptr_dn -= 2; }
735 }
736 }
737 } else {
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]);
742 if (dof.size2()) {
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); }
745 } else {
746 ptr_dn += 6;
747 }
748 if (face_id > 6) { ptr_dn -= 3; }
749 } else {
750 dof5_node_type::get_global_dof_numbering_s(ptr_dn, nodes[kk]);
751 if (dof.size2()) {
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);
754 ptr_dn++;
755 const TP d1 = dof(*ptr_dn, 0);
756 ptr_dn++;
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;
760 } else {
761 ptr_dn += 5;
762 }
763 if (face_id > 6) { ptr_dn -= 2; }
764 }
765 }
766 /* the incompatible mode are not used in the value computation
767 if (incompatible_mode && face_id < 5)
768 for (int i = 0; i != 4; ++i, ++ptr_dn) {
769 *ptr_dn = incompatible_dof_index + i;
770 inc_dof[i] = dof(*ptr_dn, 0);
771 }
772 */
773 }
774
775 // The director in the deformed configuration at each nodes.
776 TP d[nb_nodes][3];
777
778 // The derivative of the director in the deformed
779 // configuration relatively to the rotation dof for each nodes.
780 TP T[nb_nodes][3][3];
781
782 // Computation of the deformed director and this derivative at
783 // each nodes. This is done by finite rotation.
784 for (int k = 0; k != nb_nodes_f; ++k) {
785 TP norm_w;
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);
789 TP O2[6];
790 {
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;
795 O2[1] = w[0] * w[1];
796 O2[2] = w[0] * w[2];
797 O2[3] = -w2_2 - w0_2;
798 O2[4] = w[1] * w[2];
799 O2[5] = -w1_2 - w0_2;
800 }
801 TP s;
802 TP c;
803 TP cc;
804 if (ww < 0.25) {
805 TP w4 = ww * ww;
806 TP w6 = w4 * ww;
807 TP w8 = w6 * ww;
808 TP w10 = w8 * ww;
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;
812 } else {
813 s = std::sin(norm_w) / norm_w;
814 c = std::sin(norm_w * 0.5) / norm_w;
815 c = 2 * c * c;
816 cc = (norm_w - std::sin(norm_w)) / (ww * norm_w);
817 }
818 {
819 TP* const dk = d[k];
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];
827 }
828
829 if (d_value_d_dof.is_null()) { continue; }
830 s = c;
831 c = cc;
832
833 TP HT3k[3][3];
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];
844 } else {
845 {
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];
853 }
854 {
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];
862 }
863 }
864 {
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];
876 }
877 }
878 }
879
880 if (!value.is_null()) {
881 b2linalg::Vector<double> shape(nb_nodes_f);
882 if (face_id < 5) {
883 edge_shape::eval_h(&internal_coor[0], shape);
884 } else {
885 SHAPE::eval_h(&internal_coor[0], shape);
886 }
887 value.resize(3);
888 for (size_t i = 0; i != 3; ++i) {
889 TP tmp = 0;
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);
892 }
893 value[i] = tmp;
894 }
895 }
896 if (!d_value_d_icoor.is_null()) {
897 if (face_id > 4) {
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) {
903 TP tmp = 0;
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);
906 }
907 d_value_d_icoor(i, j) = tmp;
908 }
909 }
910 } else {
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) {
917 TP tmp0 = 0;
918 TP tmp1 = 0;
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]);
922 }
923 d_value_d_icoor(i, 0) = tmp0;
924 d_value_d_icoor(i, 1) = tmp1 * z_scale;
925 }
926 }
927 }
928 if (!d_value_d_dof.is_null()) {
929 dof_numbering.assign(dn, ptr_dn);
930 b2linalg::Vector<double> shape(nb_nodes_f);
931 if (face_id < 5) {
932 edge_shape::eval_h(&internal_coor[0], shape);
933 } else {
934 SHAPE::eval_h(&internal_coor[0], shape);
935 }
936 d_value_d_dof.resize(3, dof_numbering.size());
937 d_value_d_dof.set_zero();
938 size_t i_ptr = 0;
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]; }
941 i_ptr += 3;
942 if (face_id < 7) {
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];
947 }
948 }
949 }
950 }
951 if (!d_value_d_dof_dep_col.is_null() && face_id < 7) {
952 size_t i_ptr = 0;
953 if (face_id < 5) {
954 for (int k = 0; k != nb_nodes_f; ++k) {
955 if (nodes_type_6_dof[map_node[k]]) {
956 i_ptr += 3;
957 int i_max = 0;
958 for (int i = 1; i != 3; ++i) {
959 if (std::abs(R[k][i]) > std::abs(R[k][i_max])) { i_max = i; }
960 }
961 d_value_d_dof_dep_col.push_back(i_ptr + i_max);
962 i_ptr += 3;
963 } else {
964 i_ptr += 5;
965 }
966 }
967 } else {
968 for (int k = 0; k != nb_nodes_f; ++k) {
969 const int i_max = nodes_type_6_dof[map_node[k]] ? 3 : 2;
970 i_ptr += 3;
971 for (int i = 0; i != i_max; ++i, ++i_ptr) {
972 d_value_d_dof_dep_col.push_back(i_ptr);
973 }
974 }
975 }
976 }
977 }
978 }
979
980 void body_geom(
981 const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
982 b2linalg::Vector<TP>& geom, b2linalg::Matrix<TP, b2linalg::Mrectangle>& d_geom_d_icoor) {
983 if (!SHAPE_IP_OOP_SAME) { UnimplementedError() << THROW; }
984
985 // The normal director at the initial configuration for each nodes.
986 TP D[nb_nodes][3] = {};
987
988 // The coordinate at each nodes
989 TP R[nb_nodes][6] = {};
990
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]);
994 } else {
995 dof5_node_type::get_coor_s(R[k], nodes[k]);
996 }
997 }
998
999 for (int k = nb_nodes_ip; k != nb_nodes; ++k) {
1000 dof2_node_type::get_coor_s(R[k], nodes[k]);
1001 }
1002
1003 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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) { // inversed normal w.r.t. element
1009 scale_3(R[k] + 3, -1.);
1010 }
1011 std::copy(R[k] + 3, R[k] + 6, D[k]);
1012 }
1013 }
1014
1015 TP z, dz;
1016 {
1017 b2linalg::Vector<double> N(nb_nodes);
1018 SHAPE::eval_h(&internal_coor[0], N);
1019 TP e_begin, e_end;
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);
1023 }
1024
1025 if (!geom.is_null()) {
1026 geom.resize(3);
1027 b2linalg::Vector<double> shape(nb_nodes);
1028 SHAPE::eval_h(&internal_coor[0], shape);
1029 for (size_t i = 0; i != 3; ++i) {
1030 TP tmp = 0;
1031 for (int k = 0; k != nb_nodes; ++k) { tmp += shape[k] * (R[k][i] + D[k][i] * z); }
1032 geom[i] = tmp;
1033 }
1034 }
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) {
1043 TP tmp = 0;
1044 for (int k = 0; k != nb_nodes; ++k) {
1045 tmp += d_shape(j, k) * (R[k][i] + D[k][i] * z);
1046 }
1047 d_geom_d_icoor(i, j) = tmp;
1048 }
1049 {
1050 TP tmp = 0;
1051 for (int k = 0; k != nb_nodes; ++k) { tmp += shape[k] * D[k][i] * dz; }
1052 d_geom_d_icoor(i, 2) = tmp;
1053 }
1054 }
1055 }
1056 }
1057
1058 bool body_field_linear_on_dof(const std::string& field_name) override { return false; }
1059
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) {
1068 if (!SHAPE_IP_OOP_SAME) { UnimplementedError() << THROW; }
1069
1070 // The normal director at the initial configuration for each nodes.
1071 TP D[nb_nodes][3] = {};
1072
1073 // The local referential for the nodes that have 5 dofs.
1074 TP Dr[nb_nodes][2][3] = {};
1075
1076 // The coordinate at each nodes
1077 TP R[nb_nodes][6] = {};
1078
1079 int number_of_dof = incompatible_mode ? 4 : 0;
1080
1081 for (int k = 0; k != nb_nodes_ip; ++k) {
1082 if (nodes_type_6_dof[k]) {
1083 number_of_dof += 6;
1084 dof6_node_type::get_coor_s(R[k], nodes[k]);
1085 } else {
1086 number_of_dof += 5;
1087 dof5_node_type::get_coor_s(R[k], nodes[k]);
1088 }
1089 }
1090
1091 for (int k = nb_nodes_ip; k != nb_nodes; ++k) {
1092 number_of_dof += 2;
1093 dof2_node_type::get_coor_s(R[k], nodes[k]);
1094 }
1095
1096 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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) { // inversed normal w.r.t. element
1102 scale_3(R[k] + 3, -1.);
1103 }
1104 std::copy(R[k] + 3, R[k] + 6, D[k]);
1105 }
1106 }
1107
1108 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
1109 for (int k = 0; k != nb_nodes; ++k) {
1110 if (!nodes_type_6_dof[k]) {
1111 const TP e2[3] = {0, 1, 0};
1112 outer_product_3(e2, D[k], Dr[k][0]);
1113 if (normalise_3(Dr[k][0]) < 0.01) {
1114 const TP e2[3] = {1, 0, 0};
1115 outer_product_3(e2, D[k], Dr[k][0]);
1116 normalise_3(Dr[k][0]);
1117 }
1118 outer_product_3(D[k], Dr[k][0], Dr[k][1]);
1119 normalise_3(Dr[k][1]);
1120 }
1121 }
1122
1123 // The 3 translations at each nodes
1124 TP u_dof[nb_nodes][3];
1125
1126 // The 3 rotations at each nodes
1127 TP w_dof[nb_nodes][3];
1128
1129 // convert the 5 dof to 6 dof
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);
1133 } else {
1134 for (int k = 0; k != nb_nodes; ++k) {
1135 size_t dn[6];
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); }
1140 } else {
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;
1146 }
1147 }
1148 }
1149
1150 // The director in the deformed configuration at each nodes.
1151 TP d[nb_nodes][3];
1152
1153 // The derivative of the director in the deformed
1154 // configuration relatively to the rotation dof for each nodes.
1155 TP HT3[nb_nodes][3][3];
1156 TP T[nb_nodes][3][3];
1157 TP norm_w[nb_nodes];
1158
1159 // Computation of the deformed director and this derivative at
1160 // each nodes. This is done by finite rotation.
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);
1165 TP O2[6];
1166 {
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;
1176 }
1177 TP s;
1178 TP c;
1179 TP cc;
1180 if (ww < 0.25) {
1181 TP w4 = ww * ww;
1182 TP w6 = w4 * ww;
1183 TP w8 = w6 * ww;
1184 TP w10 = w8 * ww;
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;
1188 } else {
1189 s = std::sin(norm_w[k]) / norm_w[k];
1190 c = std::sin(norm_w[k] * 0.5) / norm_w[k];
1191 c = 2 * c * c;
1192 cc = (norm_w[k] - std::sin(norm_w[k])) / (ww * norm_w[k]);
1193 }
1194 {
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];
1203 }
1204
1205 if (!d_value_d_dof.is_null()) {
1206 s = c;
1207 c = cc;
1208
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];
1221 } else {
1222 {
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];
1230 }
1231 {
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];
1239 }
1240 }
1241 {
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];
1253 }
1254 }
1255 }
1256 }
1257
1258 TP z, dz;
1259 {
1260 b2linalg::Vector<double> N(nb_nodes);
1261 SHAPE::eval_h(&internal_coor[0], N);
1262 TP e_begin, e_end;
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);
1266 }
1267
1268 if (!value.is_null()) {
1269 b2linalg::Vector<double> shape(nb_nodes);
1270 SHAPE::eval_h(&internal_coor[0], shape);
1271 value.resize(3);
1272 for (size_t i = 0; i != 3; ++i) {
1273 TP tmp = 0;
1274 for (int k = 0; k != nb_nodes; ++k) {
1275 tmp += shape[k] * (u_dof[k][i] + (d[k][i] - D[k][i]) * z);
1276 }
1277 value[i] = tmp;
1278 }
1279 }
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) {
1288 TP tmp = 0;
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);
1291 }
1292 d_value_d_icoor(i, j) = tmp;
1293 }
1294 {
1295 TP tmp = 0;
1296 for (int k = 0; k != nb_nodes; ++k) {
1297 tmp += shape[k] * (d[k][i] - D[k][i]) * dz;
1298 }
1299 d_value_d_icoor(i, 2) = tmp;
1300 }
1301 }
1302 }
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();
1309 size_t i_ptr = 0;
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];
1316 }
1317 }
1318 }
1319 }
1320 }
1321
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);
1329
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);
1337
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);
1344
1345 static type_t type;
1346
1347private:
1348 static const int nb_nodes = SHAPE::num_nodes;
1349 static const int nb_nodes_ip = SHAPE_IP::num_nodes;
1350 // static const int nb_thickness_gauss = 2;
1351
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;
1357#endif
1358#else
1359 static const int nb_thickness_gauss = 2;
1360#endif
1361
1362 static const int max_nb_dof = nb_nodes * 6 + (incompatible_mode ? 4 : 0);
1363
1364 // The dof index of the incompatible dof
1365 size_t incompatible_dof_index;
1366
1367 // The nodes list.
1368 Node* nodes[nb_nodes];
1369
1370 // The property object.
1371 SolidMechanics25D<TP>* property;
1372 TP non_structural_mass;
1373
1374 bool line_load_as_moment;
1375
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);
1379 }
1380
1381 // The nodes type. If nodes_type_6_dof[i] is true the iem node is
1382 // a 6 dof node else it is a 5 dof nodes or 2 dof nodes.
1383 std::bitset<nb_nodes> nodes_type_6_dof;
1384
1385 typedef TP C_t[36];
1386 typedef SolidMechanics25D<TP>* properties_list_t[nb_gauss][nb_thickness_gauss];
1387 union {
1388 // The list of properties at each integration points for
1389 // path_dependent material.
1390 properties_list_t* properties_list;
1391
1392 // The constitutive tensors at each integration points
1393 // expressed in the contravariant reference configuration
1394 // for linear material.
1395 C_t* C_linear;
1396 };
1397
1398 struct OnePointShape {
1399 // Integration point coordinate
1400 double IntCoor[2];
1401
1402 // The shape functions at the Gauss points.
1403 double N[nb_nodes];
1404
1405 // The derivative of shape function at the Gauss points.
1406 double dN[2][nb_nodes];
1407
1408 // The shape functions at the Gauss points.
1409 double N_ip[nb_nodes_ip];
1410
1411 // The derivative of shape function at the Gauss points.
1412 double dN_ip[2][nb_nodes_ip];
1413 };
1414
1415 struct OnePoint : public OnePointShape, public MITC_TYING_POINT::TyingPointInterp {
1416 // Integration point weight
1417 double weight;
1418 };
1419
1420 // The shape functions, the derivative, the weight and the
1421 // tying point interpolation at the Gauss point.
1422 static OnePoint gauss_point[nb_gauss];
1423
1424 // The shape functions and derivative at the tying point.
1425 static OnePointShape tying_point[MITC_TYING_POINT::nb_tying_point];
1426
1427#ifdef TAYLOR_IM
1428 static OnePointShape central_point;
1429#endif
1430
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);
1434 ops.IntCoor[0] = r;
1435 ops.IntCoor[1] = s;
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); }
1441 }
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); }
1448 }
1449 }
1450
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);
1454 }
1455
1456 static void init_all_tying_point() {
1457 for (int i = 0; i != MITC_TYING_POINT::nb_tying_point; ++i) {
1458 init_shape_point(
1459 MITC_TYING_POINT::tying_point[i].r, MITC_TYING_POINT::tying_point[i].s,
1460 tying_point[i]);
1461 }
1462 }
1463
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."
1470 << THROW;
1471 }
1472 IP_ID ip_id;
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]);
1476 }
1477 }
1478
1479#ifdef CHECK_TYING_INTERPOLATION
1480 // This check is not valid for triangle element and return false error.
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")
1486 << std::endl;
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
1495 << std::endl;
1496 }
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;
1501 }
1502 }
1503 break;
1504 }
1505 }
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
1510 << std::endl;
1511 }
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;
1516 }
1517 }
1518 break;
1519 }
1520 }
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
1525 << std::endl;
1526 }
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;
1531 }
1532 }
1533 break;
1534 }
1535 }
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
1540 << std::endl;
1541 }
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;
1546 }
1547 }
1548 break;
1549 }
1550 }
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
1555 << std::endl;
1556 }
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;
1561 }
1562 break;
1563 }
1564 }
1565 }
1566 }
1567 }
1568#endif
1569
1570 struct InitElemType {
1571 InitElemType() {
1572 // we comment this because we don't when to check the code each
1573 // time b2000++ is started. Uncomment this to check the tying
1574 // point if there are modified.
1575#ifdef CHECK_TYING_INTERPOLATION
1576 ShellMITC::check_tying_interpolation();
1577#endif
1578
1579 ShellMITC::init_all_tying_point();
1580 ShellMITC::init_all_gauss_point();
1581#ifdef TAYLOR_IM
1582 init_shape_point(0, 0, central_point);
1583#endif
1584 }
1585 };
1586
1587 static const InitElemType init_elem_type;
1588};
1589
1590template <
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>
1594const int ShellMITC<
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;
1597
1598template <
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>
1602typename ShellMITC<
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
1605 ShellMITC<
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];
1608
1609template <
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>
1613typename ShellMITC<
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
1616 ShellMITC<
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];
1620
1621#ifdef TAYLOR_IM
1622template <
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>
1626typename ShellMITC<
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
1629 ShellMITC<
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;
1632#endif
1633
1634template <
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
1641 ShellMITC<
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;
1644
1645template <
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>
1649typename ShellMITC<
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
1652 ShellMITC<
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);
1660
1661} // namespace b2000
1662
1663template <
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>::
1670 get_value(
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};
1689#endif
1690#else
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};
1694#endif
1695
1696 get_dof_numbering(dof_numbering);
1697 if (!value.is_null()) {
1698 value.resize(dof_numbering.size());
1699 value.set_zero();
1700 }
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();
1706 } else {
1707 d_value_d_dof[i].resize(size_t(0), true);
1708 }
1709 }
1710 if (!d_value_d_time.is_null()) {
1711 d_value_d_time.resize(dof_numbering.size());
1712 d_value_d_time.set_zero();
1713 }
1714
1715 // The normal director at the initial configuration for each nodes.
1716 TP D[nb_nodes][3] = {};
1717 TP Dd[nb_nodes][3] = {};
1718
1719 // The local referential for the nodes that have 5 dofs.
1720 TP Dr[nb_nodes][2][3] = {};
1721
1722 // The coordinate at each node plus the normal vector or (0, 0, 0)
1723 TP R[nb_nodes][6] = {};
1724
1725 const int number_of_dof = (int)dof_numbering.size();
1726
1727 // int number_of_dof = incompatible_mode ? 4 : 0;
1728
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]);
1732 } else {
1733 dof5_node_type::get_coor_s(R[k], nodes[k]);
1734 }
1735 }
1736 for (int k = nb_nodes_ip; k != nb_nodes; ++k) { dof2_node_type::get_coor_s(R[k], nodes[k]); }
1737
1738 bool drill_stabilized_node[nb_nodes] = {};
1739
1740 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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) { // inversed normal w.r.t. element
1745 scale_3(R[k] + 3, -1.);
1746 }
1747 std::copy(R[k] + 3, R[k] + 6, D[k]);
1748 drill_stabilized_node[k] = nodes_type_6_dof[k];
1749 } else {
1750 // 6dof nodes that should not be drill-stabilized are
1751 // marked by b2add_6dof as having an empty normal vector.
1752 drill_stabilized_node[k] = false;
1753 }
1754 }
1755
1756 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
1757 for (int k = 0; k != nb_nodes; ++k) {
1758 if (!nodes_type_6_dof[k]) {
1759 const TP e2[3] = {0, 1, 0};
1760 outer_product_3(e2, D[k], Dr[k][0]);
1761 if (normalise_3(Dr[k][0]) < 0.01) {
1762 const TP e2[3] = {1, 0, 0};
1763 outer_product_3(e2, D[k], Dr[k][0]);
1764 normalise_3(Dr[k][0]);
1765 }
1766 outer_product_3(D[k], Dr[k][0], Dr[k][1]);
1767 normalise_3(Dr[k][1]);
1768 }
1769 }
1770
1771 // The 3 translations at each nodes
1772 TP u_dof[nb_nodes][3];
1773
1774 // The 3 rotations at each nodes
1775 TP w_dof[nb_nodes][3];
1776
1777 // The incompatible dof
1778 TP inc_dof[4];
1779
1780 // convert the 5 dof to 6 dof
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); }
1785 } else {
1786 size_t o = 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]);
1792 o += 6;
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];
1798 o += 5;
1799 } else {
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];
1803 o += 2;
1804 }
1805 }
1806 if (incompatible_mode) {
1807 const TP* alpha = &dof[0][dof_numbering[o]];
1808 std::copy(alpha, alpha + 4, inc_dof);
1809 }
1810 }
1811
1812 // The director in the deformed configuration at each nodes.
1813 TP d[nb_nodes][3];
1814
1815 // The derivative of the director in the deformed
1816 // configuration relatively to the rotation dof for each nodes.
1817 TP HT3[nb_nodes][3][3];
1818 TP T[nb_nodes][3][3];
1819 TP norm_w[nb_nodes];
1820
1821 // Computation of the deformed director and this derivative at
1822 // each nodes. This is done by finite rotation.
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);
1827 TP O2[6];
1828 {
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;
1838 }
1839 TP s;
1840 TP c;
1841 TP cc;
1842 if (ww < 0.25) {
1843 TP w4 = ww * ww;
1844 TP w6 = w4 * ww;
1845 TP w8 = w6 * ww;
1846 TP w10 = w8 * ww;
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;
1850 } else {
1851 s = std::sin(norm_w[k]) / norm_w[k];
1852 c = std::sin(norm_w[k] * 0.5) / norm_w[k];
1853 c = 2 * c * c;
1854 cc = (norm_w[k] - std::sin(norm_w[k])) / (ww * norm_w[k]);
1855 }
1856 {
1857 TP* const dk = d[k];
1858 const TP* const Dk = D[k];
1859 if (linear) {
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];
1867 } else {
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];
1878 }
1879 }
1880
1881 if (linear) {
1882 const TP* const dk = D[k];
1883 if (nodes_type_6_dof[k]) {
1884 T[k][0][0] = 0;
1885 T[k][0][1] = -dk[2];
1886 T[k][0][2] = dk[1];
1887 T[k][1][0] = dk[2];
1888 T[k][1][1] = 0;
1889 T[k][1][2] = -dk[0];
1890 T[k][2][0] = -dk[1];
1891 T[k][2][1] = dk[0];
1892 T[k][2][2] = 0;
1893 } else {
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];
1900 }
1901 } else {
1902 s = c;
1903 c = cc;
1904
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];
1917 } else {
1918 {
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];
1926 }
1927 {
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];
1935 }
1936 }
1937 {
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];
1949 }
1950 }
1951 }
1952 }
1953
1954 // The Green Lagrange shear strain at the tying point for the
1955 // assumed natural strain (MITC formulation).
1956 struct TyingPoint {
1957 TP ars[2][3];
1958 TP at[3];
1959 TP atrs[2][3];
1960 TP Ers0[3];
1961 TP Ers1[3];
1962 TP Erst[2];
1963 TP contravariant[mitc_gs ? 3 : 0][3];
1964 };
1965 TyingPoint tp_value[MITC_TYING_POINT::nb_tying_point];
1966
1967 // Iteration trough the tying points.
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];
1971 TP Ar[3] = {};
1972 TP As[3] = {};
1973 TP At[3] = {};
1974 TP Atr[3] = {};
1975 TP Ats[3] = {};
1976 TP ur[3] = {};
1977 TP us[3] = {};
1978 TP ut[3] = {};
1979 TP utr[3] = {};
1980 TP uts[3] = {};
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];
1991 }
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];
1995 }
1996 }
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];
2002 }
2003 }
2004 }
2005 if (linear) {
2006 for (int j = 0; j != 3; ++j) {
2007 tp.ars[0][j] = Ar[j];
2008 tp.ars[1][j] = As[j];
2009 tp.at[j] = At[j];
2010 tp.atrs[0][j] = Atr[j];
2011 tp.atrs[1][j] = Ats[j];
2012 }
2013 } else {
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];
2020 }
2021 }
2022 if (mitc_gs) {
2023 TP A[3][3];
2024 for (int j = 0; j != 3; ++j) {
2025 A[0][j] = Ar[j];
2026 A[1][j] = As[j];
2027 A[2][j] = At[j];
2028 }
2029 invert_3_3(A, tp.contravariant);
2030 }
2031 if (MITC_TYING_POINT::tying_point[l].comp_ip || mitc_gs) {
2032 // Err0
2033 tp.Ers0[0] = Ar[0] * ur[0] + Ar[1] * ur[1] + Ar[2] * ur[2];
2034 // Ess0
2035 tp.Ers0[1] = As[0] * us[0] + As[1] * us[1] + As[2] * us[2];
2036 // 2*Ers0
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];
2039 // Err1
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];
2042 // Ess1
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];
2045 // 2*Ers1
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];
2049 if (!linear) {
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];
2057 }
2058
2059 if (incompatible_mode) {
2060 // 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];
2065 if (linear) {
2066 tp.Ers0[0] += a00;
2067 tp.Ers0[1] += a11;
2068 tp.Ers0[2] += a01 + a10;
2069 } else {
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;
2073 }
2074 }
2075 }
2076 if (MITC_TYING_POINT::tying_point[l].comp_oop || mitc_gs) {
2077 // 2*Ert0
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];
2080 // 2*Est0
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];
2083 if (!linear) {
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];
2086 }
2087 }
2088 }
2089
2090 // The Green Lagrange shell strain
2091 TP E[nb_gauss][8];
2092
2093 // The PK2 stress
2094 TP S[nb_gauss][8];
2095
2096 // The derivative of the shell strain
2097 TP Bl[max_nb_dof][nb_gauss][8];
2098
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();
2111
2112 // The derivative relatively to the rotation dof of the derivative
2113 // relatively to the rotation dof of the strain time the stress.
2114 TP h[nb_nodes][3];
2115 std::fill_n(h[0], nb_nodes * 3, 0);
2116
2117 TP C_tmp[nb_gauss][36];
2118 C_t* C = 0;
2119 if (linear_mat && !linear && !property_shell_stress) {
2120 if (!C_linear) {
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;
2124 }
2125 C = C_linear;
2126 } else {
2127 if (d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]) {
2128 C = C_tmp;
2129 need_to_compute_C_linear = true;
2130 std::fill_n(C[0], 36 * nb_gauss, 0);
2131 }
2132 if (!value.is_null() || compute_stress) { std::fill_n(S[0], 8 * nb_gauss, 0); }
2133 }
2134
2135 // if (compute_value_from_dofdotdot) {
2136 // value.resize(number_of_dof);
2137 // value.set_zero();
2138 // }
2139
2140 // if (!d_value_d_dofdotdot_null) {
2141 // d_value_d_dof[2].resize(number_of_dof, true);
2142 // d_value_d_dof[2].set_zero();
2143 // }
2144
2145 // if (d_value_d_dof.size() > 1 && d_value_d_dof[1].size1() != 0) {
2146 // d_value_d_dof[1].resize(number_of_dof, true);
2147 // d_value_d_dof[1].set_zero();
2148 // }
2149
2150 TP tp_trans[mitc_gs ? nb_gauss : 0][MITC_TYING_POINT::nb_tying_point][5][5];
2151
2152 // Iteration trough the Gauss points.
2153 for (int l = 0; l != nb_gauss; ++l) {
2154 const OnePoint& N = gauss_point[l];
2155
2156 // Discretized covariant base vectors in the reference
2157 // configuration and current configuration.
2158 // A[gauss_point] = covariant base vectors
2159 TP Ar[3] = {};
2160 TP As[3] = {};
2161 TP Atr[3] = {};
2162 TP Ats[3] = {};
2163 TP ur[3] = {};
2164 TP us[3] = {};
2165 TP ut[3] = {};
2166 TP ar[3] = {}; // ar = Ar + ur
2167 TP as[3] = {}; // as = As + us
2168 TP at[3] = {}; // at = At + ut
2169 TP atr[3] = {}; // ar = Atr + utr
2170 TP ats[3] = {}; // as = Ats + uts
2171 TP utr[3] = {};
2172 TP uts[3] = {};
2173 TP At[3] = {};
2174
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];
2185 }
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];
2189 }
2190 }
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];
2196 }
2197 }
2198 }
2199
2200 if (mitc_gs) {
2201 TP A[3][3];
2202 for (int j = 0; j != 3; ++j) {
2203 A[0][j] = Ar[j];
2204 A[1][j] = As[j];
2205 A[2][j] = At[j];
2206 }
2207 for (int tk = 0; tk != MITC_TYING_POINT::nb_tying_point; ++tk) {
2208 TP trans[3][3];
2209 inner_product_3_3_NN(tp_value[tk].contravariant, A, trans);
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];
2216
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];
2223
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]);
2230
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]);
2237
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]);
2244 }
2245 }
2246
2247 // Computation of the Green Lagrange shell strain.
2248 if (gradient_container || compute_stress) {
2249 TP* const El = E[l];
2250 // Err0 and Err1
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];
2255 } else {
2256 El[0] = 0;
2257 El[3] = 0;
2258 for (int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
2259 const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
2260 if (!mitc_gs) {
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];
2263 } else {
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]);
2275 }
2276 }
2277 }
2278 // Ess0 and Ess1
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];
2283 } else {
2284 El[1] = 0;
2285 El[4] = 0;
2286 for (int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
2287 const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
2288 if (!mitc_gs) {
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];
2291 } else {
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]);
2303 }
2304 }
2305 }
2306 // 2*Ers0 and 2*Ers1
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];
2313 } else {
2314 El[2] = 0;
2315 El[5] = 0;
2316 for (int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
2317 const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
2318 if (!mitc_gs) {
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];
2321 } else {
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]);
2333 }
2334 }
2335 }
2336 // Ert0
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];
2340 } else {
2341 El[6] = 0;
2342 for (int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
2343 const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
2344 if (!mitc_gs) {
2345 El[6] += interp.weight * tp_value[interp.tp_pos].Erst[interp.e_pos];
2346 } else {
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]);
2354 }
2355 }
2356 }
2357 // Est0
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];
2361 } else {
2362 El[7] = 0;
2363 for (int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
2364 const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
2365 if (!mitc_gs) {
2366 El[7] += interp.weight * tp_value[interp.tp_pos].Erst[interp.e_pos];
2367 } else {
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]);
2375 }
2376 }
2377 }
2378 if (!linear) {
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];
2382 }
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];
2386 }
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];
2391 }
2392 if (MITC_TYING_POINT::nb_ert_comp == 0) {
2393 El[6] += ut[0] * ur[0] + ut[1] * ur[1] + ut[2] * ur[2];
2394 }
2395 if (MITC_TYING_POINT::nb_est_comp == 0) {
2396 El[7] += ut[0] * us[0] + ut[1] * us[1] + ut[2] * us[2];
2397 }
2398 }
2399 if (incompatible_mode) {
2400 // 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];
2405 if (linear) {
2406 El[0] += a00;
2407 El[1] += a11;
2408 El[2] += a01 + a10;
2409 } else {
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;
2413 }
2414 }
2415 }
2416
2417 // Calculate branch-global coordinates of the integration point.
2418 TP coor[3] = {};
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];
2423 }
2424
2425 if (gradient_container || (compute_stress && !(linear_mat && C != 0))
2426 || need_to_compute_C_linear) {
2427 if (property_shell_stress) {
2428 TP G[6][3];
2429 TP u[6][3];
2430 for (int j = 0; j != 3; ++j) {
2431 G[0][j] = Ar[j];
2432 G[1][j] = As[j];
2433 G[2][j] = At[j];
2434 G[3][j] = Atr[j];
2435 G[4][j] = Ats[j];
2436 G[5][j] = 0;
2437 u[0][j] = ur[j];
2438 u[1][j] = us[j];
2439 u[2][j] = ut[j];
2440 u[3][j] = utr[j];
2441 u[4][j] = uts[j];
2442 u[5][j] = 0;
2443 }
2444
2445 TP w = determinant_3_3(G) * gauss_point[l].weight;
2446 if (gradient_container) {
2447 // Compute and store the integration point coordinate
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,
2453 false, typeid(TP)};
2454 if (gradient_container->compute_field_value(coor_descr.name)) {
2455 gradient_container->set_field_value(coor_descr, coor);
2456 }
2457
2458 // Compute and store the displacement at the integration
2459 // point coordinate
2460 static const GradientContainer::FieldDescription disp_descr = {
2461 "DISP_IP", "dx.dy.dz", "Displacement at the point", 3, 3, 1, 3,
2462 false, typeid(TP)};
2463 if (gradient_container->compute_field_value(disp_descr.name)) {
2464 TP disp[3] = {};
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];
2469 }
2470 gradient_container->set_field_value(disp_descr, disp);
2471 }
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);
2478 }
2479 }
2480
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),
2488 coor, // bg_coordinates
2489 G,
2490 true, // use_covariant_base
2491 u, // displacement_gradient
2492 E[l], // strain
2493 (compute_stress ? S[l] : 0), // stress
2494 (need_to_compute_C_linear ? C[l] : 0), // d_stress_d...
2495 0 // d_stress_d_time
2496 );
2497
2498 if (compute_stress) {
2499 for (int i = 0; i != 8; ++i) { S[l][i] *= w; }
2500 }
2501 if (need_to_compute_C_linear) {
2502 for (int i = 0; i != 36; ++i) { C[l][i] *= w; }
2503 }
2504
2505 } else {
2506 TP e_begin;
2507 TP e_end;
2508 const int number_of_layer = property->number_of_layer();
2509 property->laminate(
2510 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin,
2511 e_end);
2512
2513 for (int layer = 0; layer != number_of_layer; ++layer) {
2514 TP l_begin;
2515 TP l_end;
2516 property->layer(
2517 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
2518 l_begin, l_end);
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) {
2523 omega = -1;
2524 } else if (std::abs(omega) < eps) {
2525 omega = 0;
2526 } else if (std::abs(omega - 1) < eps) {
2527 omega = 1;
2528 }
2529 const bool write_gradients = number_of_layer == 1 || n_omega[lt] == 0
2530 || layer == 0
2531 || (layer + 1 == number_of_layer);
2532 TP u[3][3];
2533 TP A[3][3];
2534 TP G[3][3];
2535 for (int j = 0; j != 3; ++j) {
2536 A[0][j] = Ar[j];
2537 A[1][j] = As[j];
2538 A[2][j] = At[j];
2539 u[0][j] = ur[j] + omega * utr[j];
2540 u[1][j] = us[j] + omega * uts[j];
2541 u[2][j] = ut[j];
2542 G[0][j] = Ar[j] + omega * Atr[j];
2543 G[1][j] = As[j] + omega * Ats[j];
2544 G[2][j] = A[2][j];
2545 }
2546 const TP det_G = determinant_3_3(G);
2547 assert(det_G != 0);
2548 TP* const El = E[l];
2549 const TP EE[6] = {El[0] + omega * El[3],
2550 El[1] + omega * El[4],
2551 0,
2552 El[2] + omega * El[5],
2553 El[7],
2554 El[6]};
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;
2561 // if (w < 0)
2562 // Exception()
2563 // << "The volume of the element "
2564 // << this->get_object_name() << " is negative." << THROW;
2565 if (gradient_container && write_gradients) {
2566 // Compute and store the integration point coordinate
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,
2572 false, typeid(TP)};
2573 if (gradient_container->compute_field_value(coor_descr.name)) {
2574 TP coor[3] = {};
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]);
2579 }
2580 gradient_container->set_field_value(coor_descr, coor);
2581 }
2582
2583 // Compute and store the displacement at the integration
2584 // point coordinate
2585 static const GradientContainer::FieldDescription disp_descr = {
2586 "DISP_IP", "dx.dy.dz", "Displacement at the point", 3, 3, 1, 3,
2587 false, typeid(TP)};
2588 if (gradient_container->compute_field_value(disp_descr.name)) {
2589 TP disp[3] = {};
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]));
2594 }
2595 gradient_container->set_field_value(disp_descr, disp);
2596 }
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);
2603 }
2604 }
2605
2606 TP Cc_buffer[21];
2607 TP* Cc;
2608 if (need_to_compute_C_linear) {
2609 Cc = Cc_buffer;
2610 } else {
2611 Cc = 0;
2612 }
2613 TP Ss_buffer[6];
2614 TP* Ss;
2615 if (compute_stress && !(linear_mat && C != 0)) {
2616 Ss = Ss_buffer;
2617 } else {
2618 Ss = 0;
2619 }
2620 TP density;
2621 (path_dependent_mat ? properties_list[layer][l][lt] : property)
2622 ->get_stress(
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>(
2627 nb_nodes, N.N),
2628 thickness_interpolation,
2629 coor, // bg_coordinates
2630 G /*A*/, // covariant_base
2631 true, // use_covariant_base
2632 w, // volume
2633 u, // displacement_gradient
2634 EE, // strain
2635 0, // strain_dot
2636 0, // velocity
2637 0, // acceleration
2638 Ss, // stress
2639 Cc, // d_stress_d_strain
2640 0, // d_stress_d_strain_dot
2641 0, // d_stress_d_time
2642 0, // inertia_force
2643 density);
2644
2645 if (Ss || Cc) {
2646 if (Ss) {
2647 TP* Sl = S[l];
2648 Sl[0] += w * Ss[0];
2649 Sl[1] += w * Ss[1];
2650 Sl[2] += w * Ss[3];
2651 const TP ow = omega * w;
2652 Sl[3] += ow * Ss[0];
2653 Sl[4] += ow * Ss[1];
2654 Sl[5] += ow * Ss[3];
2655 Sl[6] += w * Ss[5];
2656 Sl[7] += w * Ss[4];
2657 }
2658 if (Cc) {
2659 assert(C != 0);
2660 TP* Cl = C[l];
2661 Cl[0] += w * Cc[0];
2662 Cl[1] += w * Cc[1];
2663 Cl[2] += w * Cc[3];
2664 Cl[8] += w * Cc[6];
2665 Cl[9] += w * Cc[8];
2666 Cl[15] += w * Cc[15];
2667
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];
2678
2679 Cl[6] += w * Cc[5];
2680 Cl[7] += w * Cc[4];
2681 Cl[13] += w * Cc[10];
2682 Cl[14] += w * Cc[9];
2683 Cl[19] += w * Cc[17];
2684 Cl[20] += w * Cc[16];
2685
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];
2693
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];
2700
2701 Cl[33] += w * Cc[20];
2702 Cl[34] += w * Cc[19];
2703 Cl[35] += w * Cc[18];
2704 }
2705 }
2706 }
2707 }
2708 }
2709 }
2710
2711#ifdef TAYLOR_IM
2712 TP J_tim[incompatible_mode ? 2 : 0][incompatible_mode ? 2 : 0];
2713 if (incompatible_mode) {}
2714#endif
2715
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) {
2718 // Computation of the PK2 shell stress.
2719 blas::spmv('L', 8, 1.0, C[l], E[l], 1, 0.0, S[l], 1);
2720 }
2721
2722 if (linear) {
2723 for (int j = 0; j != 3; ++j) {
2724 ar[j] = Ar[j];
2725 as[j] = As[j];
2726 at[j] = At[j];
2727 atr[j] = Atr[j];
2728 ats[j] = Ats[j];
2729 }
2730 } else {
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];
2737 }
2738 }
2739
2740 // Compute the Bl matrix
2741 int dofi = 0;
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];
2749 } else {
2750 Blp[0] = 0;
2751 Blp[3] = 0;
2752 for (int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
2753 const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
2754 if (!mitc_gs) {
2755 if (interp.e_pos != 2) {
2756 Blp[0] +=
2757 interp.weight
2758 * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
2759 * tp_value[interp.tp_pos].ars[interp.e_pos][j];
2760 Blp[3] +=
2761 interp.weight
2762 * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
2763 * tp_value[interp.tp_pos].atrs[interp.e_pos][j];
2764 } else {
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]);
2775 }
2776 } else {
2777 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
2778 Blp[0] +=
2779 interp.weight
2780 * (tp_value[interp.tp_pos].ars[0][j]
2781 * (tpt[0]
2782 * tying_point[interp.tp_pos].dN_ip[0][k]
2783 + tpt[2]
2784 * tying_point[interp.tp_pos]
2785 .dN_ip[1][k])
2786 + tp_value[interp.tp_pos].ars[1][j]
2787 * (tpt[1]
2788 * tying_point[interp.tp_pos]
2789 .dN_ip[1][k]
2790 + tpt[2]
2791 * tying_point[interp.tp_pos]
2792 .dN_ip[0][k])
2793 + tp_value[interp.tp_pos].at[j]
2794 * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
2795 + tpt[4]
2796 * tying_point[interp.tp_pos]
2797 .dN[1][k]));
2798 Blp[3] +=
2799 interp.weight
2800 * (tp_value[interp.tp_pos].atrs[0][j]
2801 * (tpt[0]
2802 * tying_point[interp.tp_pos].dN_ip[0][k]
2803 + tpt[2]
2804 * tying_point[interp.tp_pos]
2805 .dN_ip[1][k])
2806 + tp_value[interp.tp_pos].atrs[1][j]
2807 * (tpt[1]
2808 * tying_point[interp.tp_pos]
2809 .dN_ip[1][k]
2810 + tpt[2]
2811 * tying_point[interp.tp_pos]
2812 .dN_ip[0][k]));
2813 }
2814 }
2815 }
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];
2819 } else {
2820 Blp[1] = 0;
2821 Blp[4] = 0;
2822 for (int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
2823 const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
2824 if (!mitc_gs) {
2825 if (interp.e_pos != 2) {
2826 Blp[1] +=
2827 interp.weight
2828 * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
2829 * tp_value[interp.tp_pos].ars[interp.e_pos][j];
2830 Blp[4] +=
2831 interp.weight
2832 * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
2833 * tp_value[interp.tp_pos].atrs[interp.e_pos][j];
2834 } else {
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]);
2845 }
2846 } else {
2847 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
2848 Blp[1] +=
2849 interp.weight
2850 * (tp_value[interp.tp_pos].ars[0][j]
2851 * (tpt[0]
2852 * tying_point[interp.tp_pos].dN_ip[0][k]
2853 + tpt[2]
2854 * tying_point[interp.tp_pos]
2855 .dN_ip[1][k])
2856 + tp_value[interp.tp_pos].ars[1][j]
2857 * (tpt[1]
2858 * tying_point[interp.tp_pos]
2859 .dN_ip[1][k]
2860 + tpt[2]
2861 * tying_point[interp.tp_pos]
2862 .dN_ip[0][k])
2863 + tp_value[interp.tp_pos].at[j]
2864 * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
2865 + tpt[4]
2866 * tying_point[interp.tp_pos]
2867 .dN[1][k]));
2868 Blp[4] +=
2869 interp.weight
2870 * (tp_value[interp.tp_pos].atrs[0][j]
2871 * (tpt[0]
2872 * tying_point[interp.tp_pos].dN_ip[0][k]
2873 + tpt[2]
2874 * tying_point[interp.tp_pos]
2875 .dN_ip[1][k])
2876 + tp_value[interp.tp_pos].atrs[1][j]
2877 * (tpt[1]
2878 * tying_point[interp.tp_pos]
2879 .dN_ip[1][k]
2880 + tpt[2]
2881 * tying_point[interp.tp_pos]
2882 .dN_ip[0][k]));
2883 }
2884 }
2885 }
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];
2889 } else {
2890 Blp[2] = 0;
2891 Blp[5] = 0;
2892 for (int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
2893 const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
2894 if (!mitc_gs) {
2895 if (interp.e_pos != 2) {
2896 Blp[2] +=
2897 interp.weight
2898 * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
2899 * tp_value[interp.tp_pos].ars[interp.e_pos][j];
2900 Blp[5] +=
2901 interp.weight
2902 * tying_point[interp.tp_pos].dN_ip[interp.e_pos][k]
2903 * tp_value[interp.tp_pos].atrs[interp.e_pos][j];
2904 } else {
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]);
2915 }
2916 } else {
2917 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
2918 Blp[2] +=
2919 interp.weight
2920 * (tp_value[interp.tp_pos].ars[0][j]
2921 * (tpt[0]
2922 * tying_point[interp.tp_pos].dN_ip[0][k]
2923 + tpt[2]
2924 * tying_point[interp.tp_pos]
2925 .dN_ip[1][k])
2926 + tp_value[interp.tp_pos].ars[1][j]
2927 * (tpt[1]
2928 * tying_point[interp.tp_pos]
2929 .dN_ip[1][k]
2930 + tpt[2]
2931 * tying_point[interp.tp_pos]
2932 .dN_ip[0][k])
2933 + tp_value[interp.tp_pos].at[j]
2934 * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
2935 + tpt[4]
2936 * tying_point[interp.tp_pos]
2937 .dN[1][k]));
2938 Blp[5] +=
2939 interp.weight
2940 * (tp_value[interp.tp_pos].atrs[0][j]
2941 * (tpt[0]
2942 * tying_point[interp.tp_pos].dN_ip[0][k]
2943 + tpt[2]
2944 * tying_point[interp.tp_pos]
2945 .dN_ip[1][k])
2946 + tp_value[interp.tp_pos].atrs[1][j]
2947 * (tpt[1]
2948 * tying_point[interp.tp_pos]
2949 .dN_ip[1][k]
2950 + tpt[2]
2951 * tying_point[interp.tp_pos]
2952 .dN_ip[0][k]));
2953 }
2954 }
2955 }
2956 if (MITC_TYING_POINT::nb_ert_comp == 0) {
2957 Blp[6] = N.dN[0][k] * at[j];
2958 } else {
2959 Blp[6] = 0;
2960 for (int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
2961 const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
2962 if (!mitc_gs) {
2963 Blp[6] += interp.weight
2964 * tying_point[interp.tp_pos].dN[interp.e_pos][k]
2965 * tp_value[interp.tp_pos].at[j];
2966 } else {
2967 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
2968 Blp[6] +=
2969 interp.weight
2970 * (tp_value[interp.tp_pos].ars[0][j]
2971 * (tpt[0]
2972 * tying_point[interp.tp_pos].dN_ip[0][k]
2973 + tpt[2]
2974 * tying_point[interp.tp_pos]
2975 .dN_ip[1][k])
2976 + tp_value[interp.tp_pos].ars[1][j]
2977 * (tpt[1]
2978 * tying_point[interp.tp_pos]
2979 .dN_ip[1][k]
2980 + tpt[2]
2981 * tying_point[interp.tp_pos]
2982 .dN_ip[0][k])
2983 + tp_value[interp.tp_pos].at[j]
2984 * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
2985 + tpt[4]
2986 * tying_point[interp.tp_pos]
2987 .dN[1][k]));
2988 }
2989 }
2990 }
2991 if (MITC_TYING_POINT::nb_est_comp == 0) {
2992 Blp[7] = N.dN[1][k] * at[j];
2993 } else {
2994 Blp[7] = 0;
2995 for (int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
2996 const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
2997 if (!mitc_gs) {
2998 Blp[7] += interp.weight
2999 * tying_point[interp.tp_pos].dN[interp.e_pos][k]
3000 * tp_value[interp.tp_pos].at[j];
3001 } else {
3002 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos + 3];
3003 Blp[7] +=
3004 interp.weight
3005 * (tp_value[interp.tp_pos].ars[0][j]
3006 * (tpt[0]
3007 * tying_point[interp.tp_pos].dN_ip[0][k]
3008 + tpt[2]
3009 * tying_point[interp.tp_pos]
3010 .dN_ip[1][k])
3011 + tp_value[interp.tp_pos].ars[1][j]
3012 * (tpt[1]
3013 * tying_point[interp.tp_pos]
3014 .dN_ip[1][k]
3015 + tpt[2]
3016 * tying_point[interp.tp_pos]
3017 .dN_ip[0][k])
3018 + tp_value[interp.tp_pos].at[j]
3019 * (tpt[3] * tying_point[interp.tp_pos].dN[0][k]
3020 + tpt[4]
3021 * tying_point[interp.tp_pos]
3022 .dN[1][k]));
3023 }
3024 }
3025 }
3026 }
3027 }
3028 int dofi1 = dofi;
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];
3034 } else {
3035 Blp[3] = 0;
3036 for (int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
3037 const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
3038 if (!mitc_gs) {
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];
3043 } else {
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]);
3049 }
3050 } else {
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]);
3055 Blp[3] +=
3056 interp.weight
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]
3062 + tpt[2]
3063 * tying_point[interp.tp_pos]
3064 .dN_ip[0][k]));
3065 }
3066 }
3067 }
3068 if (MITC_TYING_POINT::nb_ess_comp == 0) {
3069 Blp[4] = N.dN_ip[1][k] * as[j];
3070 } else {
3071 Blp[4] = 0;
3072 for (int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
3073 const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
3074 if (!mitc_gs) {
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];
3079 } else {
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]);
3085 }
3086 } else {
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]);
3091 Blp[4] +=
3092 interp.weight
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]
3098 + tpt[2]
3099 * tying_point[interp.tp_pos]
3100 .dN_ip[0][k]));
3101 }
3102 }
3103 }
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];
3106 } else {
3107 Blp[5] = 0;
3108 for (int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
3109 const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
3110 if (!mitc_gs) {
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];
3115 } else {
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]);
3121 }
3122 } else {
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]);
3127 Blp[5] +=
3128 interp.weight
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]
3134 + tpt[2]
3135 * tying_point[interp.tp_pos]
3136 .dN_ip[0][k]));
3137 }
3138 }
3139 }
3140 if (MITC_TYING_POINT::nb_ert_comp == 0) {
3141 Blp[6] = N.N[k] * ar[j];
3142 } else {
3143 Blp[6] = 0;
3144 for (int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
3145 const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
3146 if (!mitc_gs) {
3147 Blp[6] += interp.weight * tying_point[interp.tp_pos].N[k]
3148 * tp_value[interp.tp_pos].ars[interp.e_pos][j];
3149 } else {
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]);
3154 }
3155 }
3156 }
3157 if (MITC_TYING_POINT::nb_est_comp == 0) {
3158 Blp[7] = N.N[k] * as[j];
3159 } else {
3160 Blp[7] = 0;
3161 for (int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
3162 const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
3163 if (!mitc_gs) {
3164 Blp[7] += interp.weight * tying_point[interp.tp_pos].N[k]
3165 * tp_value[interp.tp_pos].ars[interp.e_pos][j];
3166 } else {
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]);
3171 }
3172 }
3173 }
3174 }
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];
3182 }
3183 }
3184
3185 dofi += j_max;
3186 }
3187 // incompatible mode in the Bl matrix
3188 if (incompatible_mode) {
3189 if (linear) {
3190 TP* Blp = Bl[dofi][l];
3191 Blp[0] = -2 * N.IntCoor[0];
3192 Blp[1] = 0;
3193 Blp[2] = 0;
3194 std::fill(Blp + 3, Blp + 8, 0);
3195
3196 Blp = Bl[++dofi][l];
3197 Blp[0] = 0;
3198 Blp[1] = 0;
3199 Blp[2] = -2 * N.IntCoor[1];
3200 std::fill(Blp + 3, Blp + 8, 0);
3201
3202 Blp = Bl[++dofi][l];
3203 Blp[0] = 0;
3204 Blp[1] = 0;
3205 Blp[2] = -2 * N.IntCoor[0];
3206 std::fill(Blp + 3, Blp + 8, 0);
3207
3208 Blp = Bl[++dofi][l];
3209 Blp[0] = 0;
3210 Blp[1] = -2 * N.IntCoor[1];
3211 Blp[2] = 0;
3212 std::fill(Blp + 3, Blp + 8, 0);
3213 } else {
3214 TP* Blp = Bl[dofi][l];
3215 Blp[0] = 2 * N.IntCoor[0] * (-1 + 2 * N.IntCoor[0] * inc_dof[0]);
3216 Blp[1] = 0;
3217 Blp[2] = 4 * N.IntCoor[0] * N.IntCoor[1] * inc_dof[1];
3218 std::fill(Blp + 3, Blp + 8, 0);
3219
3220 Blp = Bl[++dofi][l];
3221 Blp[0] = 0;
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);
3225
3226 Blp = Bl[++dofi][l];
3227 Blp[0] = 4 * N.IntCoor[0] * N.IntCoor[0] * inc_dof[2];
3228 Blp[1] = 0;
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);
3231
3232 Blp = Bl[++dofi][l];
3233 Blp[0] = 0;
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);
3237 }
3238 }
3239 }
3240
3241 // Append the contribution to the vector h
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];
3247 } else {
3248 for (int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
3249 const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
3250 if (!mitc_gs) {
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];
3255 } else {
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]);
3261 }
3262 } else {
3263 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
3264 h[k][j] +=
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]
3271 + tpt[2]
3272 * tying_point[interp.tp_pos].dN[0][k]));
3273 }
3274 }
3275 }
3276 if (MITC_TYING_POINT::nb_ess_comp == 0) {
3277 h[k][j] += S[l][4] * N.dN[1][k] * as[j];
3278 } else {
3279 for (int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
3280 const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
3281 if (!mitc_gs) {
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];
3286 } else {
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]);
3292 }
3293 } else {
3294 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
3295 h[k][j] +=
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]
3302 + tpt[2]
3303 * tying_point[interp.tp_pos].dN[0][k]));
3304 }
3305 }
3306 }
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]);
3309 } else {
3310 for (int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
3311 const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
3312 if (!mitc_gs) {
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];
3317 } else {
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]);
3323 }
3324 } else {
3325 const TP* tpt = tp_trans[l][interp.tp_pos][interp.e_pos];
3326 h[k][j] +=
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]
3333 + tpt[2]
3334 * tying_point[interp.tp_pos].dN[0][k]));
3335 }
3336 }
3337 }
3338 if (MITC_TYING_POINT::nb_ert_comp == 0) {
3339 h[k][j] += S[l][6] * N.N[k] * ar[j];
3340 } else {
3341 for (int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
3342 const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
3343 if (!mitc_gs) {
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];
3346 } else {
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]);
3351 }
3352 }
3353 }
3354 if (MITC_TYING_POINT::nb_est_comp == 0) {
3355 h[k][j] += S[l][7] * N.N[k] * as[j];
3356 } else {
3357 for (int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
3358 const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
3359 if (!mitc_gs) {
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];
3362 } else {
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]);
3367 }
3368 }
3369 }
3370 }
3371 }
3372 }
3373
3374 if (!d_value_d_dofdotdot_null || compute_value_from_dofdotdot) {
3375 // Compute the mass matrix
3376 TP e_begin;
3377 TP e_end;
3378 const int number_of_layer = property->number_of_layer();
3379 property->laminate(
3380 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin,
3381 e_end);
3382
3383 //
3384 // First assemble a vector of quantities for all
3385 // through-the-thickness integration points, plus the
3386 // mid-surface for the non-structural mass.
3387 //
3388 struct Data {
3389 TP omega;
3390 TP G[3][3];
3391 TP density; // resultant
3392 TP inertia_force[3];
3393 };
3394 std::vector<Data> data_vec;
3395
3396 // Iteration on each layer
3397 for (int layer = 0; layer != number_of_layer; ++layer) {
3398 TP l_begin;
3399 TP l_end;
3400 property->layer(
3401 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
3402 l_begin, l_end);
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];
3411 }
3412 data.density = property->get_density(layer) * gauss_point[l].weight
3413 * n_weight[lt] * (l_end - l_begin) * determinant_3_3(data.G);
3414 std::fill_n(data.inertia_force, 3, 0.);
3415 }
3416 }
3417 if (non_structural_mass != 0) {
3418 data_vec.push_back(Data());
3419 Data& data = data_vec.back();
3420 data.omega = 0.;
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];
3425 }
3426 data.density =
3427 gauss_point[l].weight * norm_outer_product_3(Ar, As) * non_structural_mass;
3428 std::fill_n(data.inertia_force, 3, 0.);
3429 }
3430
3431 //
3432 // Second, loop over this vector and add to the mass matrix.
3433 //
3434 for (auto data : data_vec) {
3435 TP Nl[max_nb_dof][3];
3436 std::fill_n(Nl[0], number_of_dof * 3, 0);
3437 int dofi = 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;
3442 dofi += 3;
3443 } else if (k < nb_nodes_ip) {
3444 Nl[dofi][0] = Nl[dofi + 1][1] = Nl[dofi + 2][2] = N.N_ip[k];
3445 dofi += 3;
3446 }
3447 const TP Nlk_omega = Nlk * data.omega;
3448 if (nodes_type_6_dof[k]) {
3449 Nl[dofi][0] = 0;
3450 Nl[dofi][1] = Nlk_omega * -D[k][2];
3451 Nl[dofi][2] = Nlk_omega * D[k][1];
3452
3453 Nl[++dofi][0] = Nlk_omega * D[k][2];
3454 Nl[dofi][1] = 0;
3455 Nl[dofi][2] = Nlk_omega * -D[k][0];
3456
3457 Nl[++dofi][0] = Nlk_omega * -D[k][1];
3458 Nl[dofi][1] = Nlk_omega * D[k][0];
3459 Nl[dofi][2] = 0;
3460 } else {
3461 {
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]);
3466 }
3467 {
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]);
3472 }
3473 }
3474 ++dofi;
3475 }
3476 if (incompatible_mode) {
3477 // incompatible dof
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];
3485 }
3486 }
3487 if (compute_value_from_dofdotdot) {
3488 blas::gemv(
3489 'N', 3, number_of_dof, data.density, Nl[0], 3, &dof(0, 2), 1, 1.0,
3490 data.inertia_force, 1);
3491 blas::gemv(
3492 'T', 3, number_of_dof, 1.0, Nl[0], 3, data.inertia_force, 1, 1.0,
3493 &value[0], 1);
3494 }
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);
3499 }
3500 }
3501 }
3502 }
3503 }
3504
3505 //
3506 // Calculate the artificial stiffness for the drill rotations by
3507 // integrating the in-plane shear modulus over the element volume.
3508 //
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];
3514 TP e_begin;
3515 TP e_end;
3516 property->laminate(
3517 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin,
3518 e_end);
3519
3520 // Discretized covariant base vectors in the reference
3521 // configuration.
3522 TP Ar[3] = {};
3523 TP As[3] = {};
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];
3528 }
3529 }
3530
3531 // Calculate surface associated with this integration point.
3532 TP area;
3533 {
3534 TP normal[3];
3535 outer_product_3(Ar, As, normal);
3536 area = norm_3(normal) * gauss_point[l].weight;
3537 }
3538
3539 const int number_of_layer = property->number_of_layer();
3540 for (int layer = 0; layer != number_of_layer; ++layer) {
3541 TP l_begin;
3542 TP l_end;
3543 property->layer(
3544 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
3545 l_begin, l_end);
3546 drill_stiffness_gauss[l] += property->get_inplane_shear_modulus(layer) * area
3547 * (l_end - l_begin) * drill_stiffness_factor;
3548 }
3549 }
3550 }
3551
3552 // Compute the first variation.
3553 // value = Bl' * S
3554 if (!value.is_null()) {
3555 blas::gemv(
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) {
3559 int p = 0;
3560 for (int k = 0; k != nb_nodes_ip; ++k) {
3561 p += 3;
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])
3567 * N.N[k] * N.N[k];
3568 for (int i = 0; i != 3; ++i) { value[p + i] += drill * D[k][i]; }
3569 }
3570 }
3571 p += 3;
3572 } else {
3573 p += 2;
3574 }
3575 }
3576 }
3577 }
3578
3579 // Compute the second variation.
3580 if (d_value_d_dof_flags.size() > 0 && d_value_d_dof_flags[0]) {
3581 assert(C != 0);
3582
3583 // compute d_value_d_dof= Bl' * C * Bl
3584 {
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) {
3589 {
3590 const TP* cc = C[l];
3591 for (int jj = 0; jj != 8; ++jj) {
3592 tmp[l][jj] += *cc * Bl[i][l][jj];
3593 ++cc;
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];
3597 }
3598 }
3599 }
3600 }
3601 ++i;
3602 blas::gemv(
3603 'T', nb_gauss * 8, i, 1.0, Bl[0][0], nb_gauss * 8, tmp[0], 1, 0.0, out_ptr,
3604 1);
3605 out_ptr += i;
3606 }
3607 }
3608
3609 // Add The second order derivative of strain in d_value_d_dof
3610 if (!linear) {
3611 int VD_col = 0;
3612 for (int kj = 0; kj != nb_nodes; ++kj) {
3613 const int j_max = nodes_type_6_dof[kj] ? 3 : 2;
3614 int VD_row = 0;
3615 for (int ki = 0; ki != nb_nodes; ++ki) {
3616 const int i_max = nodes_type_6_dof[ki] ? 3 : 2;
3617 TP n = 0;
3618 TP m_qrs = 0;
3619 TP m_qsr = 0;
3620 for (int l = 0; l != nb_gauss; ++l) {
3621 const OnePoint& N = gauss_point[l];
3622 TP m = 0;
3623 if (MITC_TYING_POINT::nb_err_comp == 0) {
3624 const double tmp = N.dN_ip[0][ki] * N.dN_ip[0][kj];
3625 n += S[l][0] * tmp;
3626 m += S[l][3] * tmp;
3627 } else {
3628 TP tmp = 0;
3629 for (int tk = 0; tk != MITC_TYING_POINT::nb_err_comp; ++tk) {
3630 const typename MITC_TYING_POINT::Interp& interp = N.err[tk];
3631 if (!mitc_gs) {
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];
3636 } else {
3637 tmp +=
3638 interp.weight
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]);
3643 }
3644 } else {
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]
3648 * (tpt[0]
3649 * tying_point[interp.tp_pos]
3650 .dN_ip[0][kj]
3651 + tpt[2]
3652 * tying_point[interp.tp_pos]
3653 .dN_ip[1][kj])
3654 + tying_point[interp.tp_pos].dN_ip[0][kj]
3655 * (tpt[1]
3656 * tying_point[interp.tp_pos]
3657 .dN_ip[0][ki]
3658 + tpt[2]
3659 * tying_point[interp.tp_pos]
3660 .dN_ip[1][ki]));
3661 }
3662 }
3663 n += S[l][0] * tmp;
3664 m += S[l][3] * tmp;
3665 }
3666 if (MITC_TYING_POINT::nb_ess_comp == 0) {
3667 const double tmp = N.dN_ip[1][ki] * N.dN_ip[1][kj];
3668 n += S[l][1] * tmp;
3669 m += S[l][4] * tmp;
3670 } else {
3671 TP tmp = 0;
3672 for (int tk = 0; tk != MITC_TYING_POINT::nb_ess_comp; ++tk) {
3673 const typename MITC_TYING_POINT::Interp& interp = N.ess[tk];
3674 if (!mitc_gs) {
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];
3679 } else {
3680 tmp +=
3681 interp.weight
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]);
3686 }
3687 } else {
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]
3691 * (tpt[0]
3692 * tying_point[interp.tp_pos]
3693 .dN_ip[0][kj]
3694 + tpt[2]
3695 * tying_point[interp.tp_pos]
3696 .dN_ip[1][kj])
3697 + tying_point[interp.tp_pos].dN_ip[0][kj]
3698 * (tpt[1]
3699 * tying_point[interp.tp_pos]
3700 .dN_ip[0][ki]
3701 + tpt[2]
3702 * tying_point[interp.tp_pos]
3703 .dN_ip[1][ki]));
3704 }
3705 }
3706 n += S[l][1] * tmp;
3707 m += S[l][4] * tmp;
3708 }
3709 if (MITC_TYING_POINT::nb_ers_comp == 0) {
3710 const double tmp =
3711 N.dN_ip[0][ki] * N.dN_ip[1][kj] + N.dN_ip[1][ki] * N.dN_ip[0][kj];
3712 n += S[l][2] * tmp;
3713 m += S[l][5] * tmp;
3714 } else {
3715 TP tmp = 0;
3716 for (int tk = 0; tk != MITC_TYING_POINT::nb_ers_comp; ++tk) {
3717 const typename MITC_TYING_POINT::Interp& interp = N.ers[tk];
3718 if (!mitc_gs) {
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];
3723 } else {
3724 tmp +=
3725 interp.weight
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]);
3730 }
3731 } else {
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]
3735 * (tpt[0]
3736 * tying_point[interp.tp_pos]
3737 .dN_ip[0][kj]
3738 + tpt[2]
3739 * tying_point[interp.tp_pos]
3740 .dN_ip[1][kj])
3741 + tying_point[interp.tp_pos].dN_ip[0][kj]
3742 * (tpt[1]
3743 * tying_point[interp.tp_pos]
3744 .dN_ip[0][ki]
3745 + tpt[2]
3746 * tying_point[interp.tp_pos]
3747 .dN_ip[1][ki]));
3748 }
3749 }
3750 n += S[l][2] * tmp;
3751 m += S[l][5] * tmp;
3752 }
3753 TP qrs = 0;
3754 TP qsr = 0;
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];
3758 } else {
3759 for (int tk = 0; tk != MITC_TYING_POINT::nb_ert_comp; ++tk) {
3760 const typename MITC_TYING_POINT::Interp& interp = N.ert[tk];
3761 if (!mitc_gs) {
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];
3768 } else {
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]);
3778 }
3779 }
3780 }
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];
3784 } else {
3785 for (int tk = 0; tk != MITC_TYING_POINT::nb_est_comp; ++tk) {
3786 const typename MITC_TYING_POINT::Interp& interp = N.est[tk];
3787 if (!mitc_gs) {
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];
3794 } else {
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]);
3804 }
3805 }
3806 }
3807 m_qrs += m + qrs;
3808 m_qsr += m + qsr;
3809 }
3810 m_qrs *= 0.5;
3811 m_qsr *= 0.5;
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];
3815 }
3816 }
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];
3820 }
3821 }
3822 if (ki >= kj) {
3823 for (int i = 0; i != 3; ++i) {
3824 d_value_d_dof[0](VD_row + i, VD_col + i) += n;
3825 }
3826 }
3827 if (ki == kj) {
3828 TP bki[3];
3829 outer_product_3(d[ki], h[ki], bki);
3830 TP bki_wki =
3831 bki[0] * w_dof[ki][0] + bki[1] * w_dof[ki][1] + bki[2] * w_dof[ki][2];
3832 TP c10;
3833 TP c11;
3834 TP c3;
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;
3840 } else {
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);
3849 }
3850
3851 c10 = c10 * bki_wki
3852 - (d[ki][0] * h[ki][0] + d[ki][1] * h[ki][1] + d[ki][2] * h[ki][2]);
3853 c11 *= bki_wki;
3854 for (int i = 0; i != 3; ++i) { bki[i] = -c3 * bki[i] + c11 * w_dof[ki][i]; }
3855 TP M[3][3];
3856 for (int i = 0; i != 3; ++i) {
3857 for (int j = 0; j != 3; ++j) {
3858 M[j][i] = .5
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]);
3861 }
3862 M[i][i] += c10;
3863 }
3864 for (int j = 0; j != i_max; ++j) {
3865 TP MHT3[3];
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];
3869 }
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];
3874 }
3875 }
3876 }
3877 VD_row += 3 + i_max;
3878 }
3879 VD_col += 3 + j_max;
3880 }
3881 if (incompatible_mode) {
3882 // incompatible mode
3883 TP a00 = 0;
3884 TP a01 = 0;
3885 TP a11 = 0;
3886 TP a22 = 0;
3887 TP a23 = 0;
3888 TP a33 = 0;
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];
3900 }
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;
3907 }
3908 }
3909
3910 if (drill_stiffness_factor != 0) {
3911 int p = 0;
3912 for (int k = 0; k != nb_nodes_ip; ++k) {
3913 p += 3;
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]
3922 * N.N[k];
3923 }
3924 }
3925 }
3926 }
3927 p += 3;
3928 } else {
3929 p += 2;
3930 }
3931 }
3932 }
3933 }
3934}
3935
3936template <
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) {
3950 if (!dofdot.is_null()) { UnimplementedError() << THROW; }
3951
3952 if (!dofdotdot.is_null()) { UnimplementedError() << THROW; }
3953
3954 if (edge < 0 || edge >= SHAPE::num_edge) {
3955 Exception() << "Invalid edge number " << edge + 1 << " for element "
3956 << this->get_object_name() << "." << THROW;
3957 }
3958
3959 if (!d_discretised_field_d_dof.is_null()) { UnimplementedError() << THROW; }
3960
3961 if (discretised_field.is_null()) { return; }
3962
3963 typedef typename SHAPE::edge_shape_type_0 edge_shape;
3964
3965 // Get the dof-numbering for the displacement dofs of the nodes
3966 // along the edge.
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);
3971 size_t ndn[6];
3972 if (nodes_type_6_dof[ii]) {
3973 dof6_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
3974 } else {
3975 dof5_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
3976 }
3977 std::copy(ndn, ndn + 3, edn + 3 * i);
3978 }
3979
3980 // Get the dof-numbering for the all dofs of the nodes along the
3981 // edge.
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);
3987 size_t ndn[6];
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);
3992 } else {
3993 dof5_node_type::get_global_dof_numbering_s(ndn, nodes[ii]);
3994 adn.insert(adn.end(), ndn, ndn + 5);
3995 }
3996 }
3997 dpos.push_back(adn.size());
3998 assert(dpos.size() == edge_shape::num_nodes + 1);
3999 assert(dpos.front() == 0 && dpos.back() == adn.size());
4000
4001 discretised_field.resize(adn.size());
4002 discretised_field.set_zero();
4003 if (!dof_numbering.is_null()) { dof_numbering = adn; }
4004
4005 // Extract the element-local, undeformed, and deformed coordinates
4006 // of the nodes along the edge.
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]);
4017 } else {
4018 dof5_node_type::get_coor_s(enode_x_0[i], nodes[ii]);
4019 }
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]]; }
4023 }
4024 }
4025
4026 // Loop over all integration points along the edge.
4027 IS_L ischeme;
4028 ischeme.set_order(edge_shape::order);
4029 for (int l = 0; l != ischeme.get_num(); ++l) {
4030 double r;
4031 double weight;
4032 ischeme.get_point(l, r, weight);
4033
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);
4038
4039 // Get the edge director.
4040 TP dr_0[3] = {};
4041 TP dr_d[3] = {};
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];
4046 }
4047 }
4048
4049 // Compute the length associated with this integration point.
4050 const TP length = normalise_3(dr_0) * weight;
4051 normalise_3(dr_d);
4052
4053 // Use the edge shape functions to interpolate the
4054 // coordinates of this integration point.
4055 double xi[3] = {};
4056 TP x_0[3] = {};
4057 TP x_d[3] = {};
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];
4063 }
4064 }
4065
4066 // Obtain the value from the Field object.
4067 b2linalg::Vector<TP, b2linalg::Vdense> value;
4068 field.get_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);
4074 if (value.size() != 3) { UnimplementedError() << THROW; }
4075
4076 // Distribute value to nodes (discretize).
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];
4081 }
4082 }
4083 } else {
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];
4087 }
4088 }
4089 }
4090 }
4091}
4092
4093template <
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) {
4107 if (!SHAPE_IP_OOP_SAME) { UnimplementedError() << THROW; }
4108
4109 if (!dofdot.is_null()) { UnimplementedError() << THROW; }
4110
4111 if (!dofdotdot.is_null()) { UnimplementedError() << THROW; }
4112
4113 if (face < 4) {
4114 UnimplementedError() << "MITC shell element " << this->get_object_name()
4115 << " cannot integrate face/edge fields over edges" << THROW;
4116 }
4117
4118 if (!d_discretised_field_d_dof.is_null()) { UnimplementedError() << THROW; }
4119
4120 // The normal director at the initial configuration for each nodes.
4121 TP D[nb_nodes][3] = {};
4122 TP Dd[nb_nodes][3] = {};
4123
4124 // The local referential for the nodes that have 5 dofs.
4125 TP Dr[nb_nodes][2][3] = {};
4126
4127 // The coordinate at each nodes
4128 TP R[nb_nodes][6] = {};
4129
4130 int number_of_dof = incompatible_mode ? 4 : 0;
4131
4132 for (int k = 0; k != nb_nodes_ip; ++k) {
4133 if (nodes_type_6_dof[k]) {
4134 number_of_dof += 6;
4135 dof6_node_type::get_coor_s(R[k], nodes[k]);
4136 } else {
4137 number_of_dof += 5;
4138 dof5_node_type::get_coor_s(R[k], nodes[k]);
4139 }
4140 }
4141
4142 for (int k = nb_nodes_ip; k != nb_nodes; ++k) {
4143 number_of_dof += 2;
4144 dof2_node_type::get_coor_s(R[k], nodes[k]);
4145 }
4146
4147 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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) { // inversed normal w.r.t. element
4153 scale_3(R[k] + 3, -1.);
4154 }
4155 std::copy(R[k] + 3, R[k] + 6, D[k]);
4156 }
4157 }
4158
4159 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
4160 for (int k = 0; k != nb_nodes; ++k) {
4161 if (!nodes_type_6_dof[k]) {
4162 const TP e2[3] = {0, 1, 0};
4163 outer_product_3(e2, D[k], Dr[k][0]);
4164 if (normalise_3(Dr[k][0]) < 0.01) {
4165 const TP e2[3] = {1, 0, 0};
4166 outer_product_3(e2, D[k], Dr[k][0]);
4167 normalise_3(Dr[k][0]);
4168 }
4169 outer_product_3(D[k], Dr[k][0], Dr[k][1]);
4170 normalise_3(Dr[k][1]);
4171 }
4172 }
4173
4174 // The 3 translations at each nodes
4175 TP u_dof[nb_nodes][3];
4176
4177 // The 3 rotations at each nodes
4178 TP w_dof[nb_nodes][3];
4179
4180 // The incompatible dof
4181 TP inc_dof[4];
4182
4183 // convert the 5 dof to 6 dof
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); }
4188 } else {
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]);
4193 alpha += 3;
4194 std::copy(alpha, alpha + 3, w_dof[k]);
4195 alpha += 3;
4196 } else {
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];
4201 alpha += 5;
4202 }
4203 }
4204 }
4205
4206 // The director in the deformed configuration at each nodes.
4207 TP d[nb_nodes][3];
4208
4209 // The derivative of the director in the deformed
4210 // configuration relatively to the rotation dof for each nodes.
4211 TP norm_w[nb_nodes];
4212
4213 // Computation of the deformed director and this derivative at
4214 // each nodes. This is done by finite rotation.
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);
4219 TP O2[6];
4220 {
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;
4230 }
4231 TP s;
4232 TP c;
4233 if (ww < 0.25) {
4234 TP w4 = ww * ww;
4235 TP w6 = w4 * ww;
4236 TP w8 = w6 * ww;
4237 TP w10 = w8 * ww;
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;
4240 } else {
4241 s = std::sin(norm_w[k]) / norm_w[k];
4242 c = std::sin(norm_w[k] * 0.5) / norm_w[k];
4243 c = 2 * c * c;
4244 }
4245 {
4246 TP* const dk = d[k];
4247 const TP* const Dk = D[k];
4248 if (0 /*linear*/) {
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];
4256 } else {
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];
4267 }
4268 }
4269 }
4270
4271 if (!discretised_field.is_null()) {
4272 discretised_field.resize(number_of_dof);
4273 discretised_field.set_zero();
4274 }
4275
4276 if (!dof_numbering.is_null()) { get_dof_numbering(dof_numbering); }
4277
4278 // Iteration over the Gauss points.
4279 for (int l = 0; l != nb_gauss; ++l) {
4280 const OnePoint& N = gauss_point[l];
4281 TP e_begin;
4282 TP e_end;
4283 property->laminate(
4284 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin, e_end);
4285
4286 // Discretized covariant base vectors in the reference
4287 // configuration and current configuration.
4288 TP Ar[3] = {};
4289 TP As[3] = {};
4290 TP At[3] = {};
4291 TP Atr[3] = {};
4292 TP Ats[3] = {};
4293 TP ur[3] = {};
4294 TP us[3] = {};
4295 TP ut[3] = {};
4296 TP utr[3] = {};
4297 TP uts[3] = {};
4298
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];
4311 }
4312 }
4313
4314 const double omega = (face == 4 ? -0.5 : (face == 5 ? +0.5 : 0.0));
4315
4316 // Compute surface "directors" and the area of the
4317 // undeformed surface.
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];
4325 }
4326 TP area;
4327 {
4328 TP normal[3];
4329 outer_product_3(d_coor_d_lcoor[0], d_coor_d_lcoor[1], normal);
4330 area = norm_3(normal) * gauss_point[l].weight;
4331 }
4332 normalise_3(d_coor_d_lcoor[0]);
4333 normalise_3(d_coor_d_lcoor[1]);
4334 normalise_3(d_coor_deformed_d_lcoor[0]);
4335 normalise_3(d_coor_deformed_d_lcoor[1]);
4336
4337 // Element-local coordinates of integration point.
4338 const double xi[3] = {N.IntCoor[0], N.IntCoor[1], omega / (e_end - e_begin) * 2};
4339
4340 // Interpolate coordinates of initial configuration.
4341 TP x_0[3] = {};
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]);
4346 }
4347
4348 // Interpolate coordinates of deformed configuration.
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]));
4354 }
4355
4356 // Obtain the value from the Field object.
4357 b2linalg::Vector<TP, b2linalg::Vdense> value;
4358 field.get_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]),
4365 value);
4366 if (value.size() != 3) { UnimplementedError() << THROW; }
4367
4368 // Scale load.
4369 if (0) {
4370 TP x = x_0[0];
4371 TP y = x_0[1];
4372 TP z = x_0[2];
4373 if (1) {
4374 // non-inhibited cylindrical shell, periodic loading
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;
4380 value *= scale;
4381 } else {
4382 const TP pi = std::acos(TP(-1));
4383 if (0) { // skewed
4384 const TP ca = 0.8660254037844387;
4385 const TP sa = 0.5;
4386 y = y / ca;
4387 x = x - y * sa;
4388 }
4389 const TP scale = std::sin(pi * x) * std::sin(pi * y);
4390 value *= scale;
4391 }
4392 }
4393
4394 TP Nl[max_nb_dof][3];
4395 std::fill_n(Nl[0], number_of_dof * 3, 0);
4396 int dofi = 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;
4400 dofi += 3;
4401 const TP Nlk_omega = Nlk * omega;
4402 if (nodes_type_6_dof[k]) {
4403 Nl[dofi][0] = 0;
4404 Nl[dofi][1] = Nlk_omega * -D[k][2];
4405 Nl[dofi][2] = Nlk_omega * D[k][1];
4406
4407 Nl[++dofi][0] = Nlk_omega * D[k][2];
4408 Nl[dofi][1] = 0;
4409 Nl[dofi][2] = Nlk_omega * -D[k][0];
4410
4411 Nl[++dofi][0] = Nlk_omega * -D[k][1];
4412 Nl[dofi][1] = Nlk_omega * D[k][0];
4413 Nl[dofi][2] = 0;
4414 } else {
4415 {
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]);
4420 }
4421 {
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]);
4426 }
4427 }
4428 ++dofi;
4429 }
4430
4431 blas::gemv(
4432 'T', 3, number_of_dof, area, Nl[0], 3, &value[0], 1, 1.0, &discretised_field[0], 1);
4433 }
4434}
4435
4436template <
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) {
4450 if (!SHAPE_IP_OOP_SAME) { UnimplementedError() << THROW; }
4451
4452 if (!dofdot.is_null()) { UnimplementedError() << THROW; }
4453
4454 if (!dofdotdot.is_null()) { UnimplementedError() << THROW; }
4455
4456 if (!d_discretised_field_d_dof.is_null()) { UnimplementedError() << THROW; }
4457
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};
4468#endif
4469#else
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};
4473#endif
4474
4475 // The normal director at the initial configuration for each nodes.
4476 TP D[nb_nodes][3] = {};
4477 TP Dd[nb_nodes][3] = {};
4478
4479 // The local referential for the nodes that have 5 dofs.
4480 TP Dr[nb_nodes][2][3] = {};
4481
4482 // The coordinate at each nodes
4483 TP R[nb_nodes][6] = {};
4484
4485 int number_of_dof = incompatible_mode ? 4 : 0;
4486
4487 for (int k = 0; k != nb_nodes_ip; ++k) {
4488 if (nodes_type_6_dof[k]) {
4489 number_of_dof += 6;
4490 dof6_node_type::get_coor_s(R[k], nodes[k]);
4491 } else {
4492 number_of_dof += 5;
4493 dof5_node_type::get_coor_s(R[k], nodes[k]);
4494 }
4495 }
4496
4497 for (int k = nb_nodes_ip; k != nb_nodes; ++k) {
4498 number_of_dof += 2;
4499 dof2_node_type::get_coor_s(R[k], nodes[k]);
4500 }
4501
4502 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
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) { // inversed normal w.r.t. element
4508 scale_3(R[k] + 3, -1.);
4509 }
4510 std::copy(R[k] + 3, R[k] + 6, D[k]);
4511 }
4512 }
4513
4514 // Modification here must also be reported in the ShellMITC_Volume_Coupling element.
4515 for (int k = 0; k != nb_nodes; ++k) {
4516 if (!nodes_type_6_dof[k]) {
4517 const TP e2[3] = {0, 1, 0};
4518 outer_product_3(e2, D[k], Dr[k][0]);
4519 if (normalise_3(Dr[k][0]) < 0.01) {
4520 const TP e2[3] = {1, 0, 0};
4521 outer_product_3(e2, D[k], Dr[k][0]);
4522 normalise_3(Dr[k][0]);
4523 }
4524 outer_product_3(D[k], Dr[k][0], Dr[k][1]);
4525 normalise_3(Dr[k][1]);
4526 }
4527 }
4528
4529 // The 3 translations at each nodes
4530 TP u_dof[nb_nodes][3];
4531
4532 // The 3 rotations at each nodes
4533 TP w_dof[nb_nodes][3];
4534
4535 // The incompatible dof
4536 TP inc_dof[4];
4537
4538 // convert the 5 dof to 6 dof
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); }
4543 } else {
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]);
4548 alpha += 3;
4549 std::copy(alpha, alpha + 3, w_dof[k]);
4550 alpha += 3;
4551 } else {
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];
4556 alpha += 5;
4557 }
4558 }
4559 }
4560
4561 // The director in the deformed configuration at each nodes.
4562 TP d[nb_nodes][3];
4563
4564 // The derivative of the director in the deformed
4565 // configuration relatively to the rotation dof for each nodes.
4566 TP norm_w[nb_nodes];
4567
4568 // Computation of the deformed director and this derivative at
4569 // each nodes. This is done by finite rotation.
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);
4574 TP O2[6];
4575 {
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;
4585 }
4586 TP s;
4587 TP c;
4588 if (ww < 0.25) {
4589 TP w4 = ww * ww;
4590 TP w6 = w4 * ww;
4591 TP w8 = w6 * ww;
4592 TP w10 = w8 * ww;
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;
4595 } else {
4596 s = std::sin(norm_w[k]) / norm_w[k];
4597 c = std::sin(norm_w[k] * 0.5) / norm_w[k];
4598 c = 2 * c * c;
4599 }
4600 {
4601 TP* const dk = d[k];
4602 const TP* const Dk = D[k];
4603 if (0 /*linear*/) { // needs API change
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];
4611 } else {
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];
4622 }
4623 }
4624 }
4625
4626 if (!discretised_field.is_null()) {
4627 discretised_field.resize(number_of_dof);
4628 discretised_field.set_zero();
4629 }
4630
4631 if (!dof_numbering.is_null()) { get_dof_numbering(dof_numbering); }
4632
4633 // Iteration over the Gauss points.
4634 for (int l = 0; l != nb_gauss; ++l) {
4635 const OnePoint& N = gauss_point[l];
4636 TP e_begin;
4637 TP e_end;
4638 const int number_of_layer = property->number_of_layer();
4639 property->laminate(
4640 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), e_begin, e_end);
4641
4642 // Discretized covariant base vectors in the reference
4643 // configuration and current configuration.
4644 TP Ar[3] = {};
4645 TP As[3] = {};
4646 TP At[3] = {};
4647 TP Atr[3] = {};
4648 TP Ats[3] = {};
4649 TP ur[3] = {};
4650 TP us[3] = {};
4651 TP ut[3] = {};
4652 TP utr[3] = {};
4653 TP uts[3] = {};
4654
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];
4667 }
4668 }
4669
4670 for (int layer = 0; layer != number_of_layer; ++layer) {
4671 TP l_begin;
4672 TP l_end;
4673 property->layer(
4674 b2linalg::Vector<double, b2linalg::Vdense_constref>(nb_nodes, N.N), layer,
4675 l_begin, l_end);
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);
4678 TP G_0[3][3];
4679 TP G_d[3][3];
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];
4683 G_0[2][j] = At[j];
4684 G_d[0][j] = ur[j] + omega * utr[j];
4685 G_d[1][j] = us[j] + omega * uts[j];
4686 G_d[2][j] = ut[j];
4687 }
4688
4689 // Element-local coordinates of integration point.
4690 const double xi[3] = {N.IntCoor[0], N.IntCoor[1], omega / (e_end - e_begin) * 2};
4691
4692 // Interpolate coordinates of initial configuration.
4693 TP x_0[3] = {};
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]);
4698 }
4699
4700 // Interpolate coordinates of deformed configuration.
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]));
4706 }
4707
4708 // Obtain the value from the Field object.
4709 b2linalg::Vector<TP, b2linalg::Vdense> value;
4710 field.get_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);
4716 if (value.size() != 3) { UnimplementedError() << THROW; }
4717 TP mult = 1.;
4718 TP add = 0.;
4719 if (field.get_multiply_with_density()) {
4720 mult = property->get_density(layer);
4721 add = (layer == 0 && lt == 1 ? non_structural_mass : TP(0));
4722 }
4723
4724 // Compute the volume associated with this integration point.
4725 const TP det_G = determinant_3_3(G_0);
4726 const TP volume = n_weight[lt] * (l_end - l_begin) * det_G * gauss_point[l].weight;
4727 const TP area = norm_outer_product_3(Ar, As) * gauss_point[l].weight;
4728
4729 TP Nl[max_nb_dof][3];
4730 std::fill_n(Nl[0], number_of_dof * 3, 0);
4731 int dofi = 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;
4735 dofi += 3;
4736 const TP Nlk_omega = Nlk * omega;
4737 if (nodes_type_6_dof[k]) {
4738 Nl[dofi][0] = 0;
4739 Nl[dofi][1] = Nlk_omega * -D[k][2];
4740 Nl[dofi][2] = Nlk_omega * D[k][1];
4741
4742 Nl[++dofi][0] = Nlk_omega * D[k][2];
4743 Nl[dofi][1] = 0;
4744 Nl[dofi][2] = Nlk_omega * -D[k][0];
4745
4746 Nl[++dofi][0] = Nlk_omega * -D[k][1];
4747 Nl[dofi][1] = Nlk_omega * D[k][0];
4748 Nl[dofi][2] = 0;
4749 } else {
4750 {
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]);
4755 }
4756 {
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]);
4761 }
4762 }
4763 ++dofi;
4764 }
4765
4766 blas::gemv(
4767 'T', 3, number_of_dof, (volume * mult + area * add), Nl[0], 3, &value[0], 1,
4768 1.0, &discretised_field[0], 1);
4769 }
4770 }
4771 }
4772}
4773
4774namespace b2000 {
4775// element to couple the first face of a MITC shell element to a
4776// node of a volume element
4777template <
4778 typename TP, typename SHAPE, typename DIRECTOR_ESTIMATION, typename EDGE, bool FREEZ = false>
4779class ShellMITC_Volume_Coupling : public b2000::TypedElement<TP> {
4780public:
4781 using dof5_node_type = node::HNode<
4782 coordof::Coor<
4783 coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
4784 coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
4785 coordof::Dof<
4786 coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
4787 coordof::DTrans<coordof::RX, coordof::RY>>>;
4788
4789 using dof6_node_type = node::HNode<
4790 coordof::Coor<
4791 coordof::Trans<coordof::X, coordof::Y, coordof::Z>,
4792 coordof::Trans<coordof::NX, coordof::NY, coordof::NZ>>,
4793 coordof::Dof<
4794 coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>,
4795 coordof::DTrans<coordof::RX, coordof::RY, coordof::RZ>>>;
4796
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>>>;
4800
4801 using type_t = ObjectTypeComplete<ShellMITC_Volume_Coupling, typename TypedElement<TP>::type_t>;
4802
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;
4806 }
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])) {
4810 Exception() << "MITC shell element " << this->get_object_name()
4811 << " cannot accept node " << nodes[i]->get_object_name()
4812 << " because it does not have the dofs UX,UY,UZ,RX,RY." << THROW;
4813 }
4814 nodes_type_6_dof[i] = dof6_node_type::is_dof_subset_s(nodes[i]);
4815 }
4816 nodes[nb_nodes] = nodes_.second[nb_nodes];
4817 if (!volume_node_type::is_dof_subset_s(nodes[nb_nodes])) {
4818 Exception() << "MITC shell element " << this->get_object_name()
4819 << " cannot accept as last node " << nodes[nb_nodes]->get_object_name()
4820 << " because it does "
4821 "not have the dofs UX,UY,UZ."
4822 << THROW;
4823 }
4824 }
4825
4826 std::pair<int, Node* const*> get_nodes() const {
4827 return std::pair<int, Node* const*>(nb_nodes + 1, nodes);
4828 }
4829
4830 void get_dof_numbering(b2linalg::Index& dof_numbering) {
4831 int s = 0;
4832 for (int k = 0; k != EDGE::nb_nodes_of; ++k) {
4833 s += nodes_type_6_dof[EDGE::nodes_of_id[k]] ? 6 : 5;
4834 }
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]]);
4840 } else {
4841 ptr = dof5_node_type::get_global_dof_numbering_s(ptr, nodes[EDGE::nodes_of_id[k]]);
4842 }
4843 }
4844 ptr = volume_node_type::get_global_dof_numbering_s(ptr, nodes[nb_nodes]);
4845 }
4846
4847 void set_property(ElementProperty* property_) { property = property_; }
4848
4849 const ElementProperty* get_property() const { return property; }
4850
4851 void init(Model& model) {}
4852
4853 // void get_value(
4854 // const b2linalg::Matrix<TP, b2linalg::Mrectangle_constref>& dof,
4855 // const TP time,
4856 // const TP delta_time,
4857 // const EquilibriumSolution equilibrium_solution,
4858 // const bool linear,
4859 // Model& model,
4860 // b2linalg::Index& dof_numbering,
4861 // b2linalg::Vector<TP, b2linalg::Vdense>& value,
4862 // const std::vector<bool>& d_value_d_dof_flags,
4863 // std::vector<b2linalg::Matrix<TP, b2linalg::Mpacked> >& d_value_d_dof,
4864 // b2linalg::Vector<TP, b2linalg::Vdense>& d_value_d_time,
4865 // GradientContainer* gradient_container,
4866 // SolverHints* solver_hints) {
4867 // if (!value.is_null()) {
4868 // size_t number_of_dof = 3;
4869 // for (int k = 0; k != EDGE::nb_nodes_of; ++k)
4870 // if (nodes_type_6_dof[EDGE::nodes_of_id[k]])
4871 // number_of_dof += 6;
4872 // else
4873 // number_of_dof += 5;
4874 // value.resize(number_of_dof);
4875 // value.set_zero();
4876 // }
4877 // if (d_value_d_dof.size() != 0 && d_value_d_dof[0].size1() != 0)
4878 // d_value_d_dof[0].set_zero();
4879 // }
4880
4881 std::pair<int, Element::VariableInfo> get_constraint_info() {
4882 return std::pair<int, Element::VariableInfo>(FREEZ ? 2 : 3, Element::nonlinear);
4883 }
4884
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);
4891
4892 static type_t type;
4893
4894private:
4895 static const int nb_nodes = SHAPE::num_nodes;
4896 ElementProperty* property;
4897 Node* nodes[nb_nodes + 1];
4898
4899 // The nodes type. If nodes_type_6_dof[i] is true the iem node is
4900 // a 6 dof node else it is a 5 dof nodes or 2 dof nodes.
4901 std::bitset<nb_nodes> nodes_type_6_dof;
4902
4903 // compute the point in the quadratic edge defined by a, b and c that
4904 // minimise the distance between d and the quadratic edge and return
4905 // the weight associated to a, b and c to interpolate this point.
4906 // a---------c-----b
4907 // d
4908 void get_p_weight(
4909 const TP a_coor[3], const TP b_coor[3], const TP c_coor[3], const TP d_coor[3], TP& va,
4910 TP& vb, TP& vc) {
4911 TP coef[4] = {};
4912 int i, j;
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];
4926 }
4927
4928 TP x[3];
4929 int nsol;
4930 solve_cubic_equation(coef[0], coef[1], coef[2], coef[3], nsol, x[0], x[1], x[2]);
4931 TP r = 0;
4932 if (nsol == 1) {
4933 r = x[0];
4934 if (r < -1) { r = -1; }
4935 if (r > 1) { r = 1; }
4936 } else {
4937 TP lr2 = 1e100;
4938 TP l_border = 1e100;
4939 for (j = 0; j != nsol; ++j) {
4940 TP lr2_tmp = 0;
4941 TP l_border_tmp = 0;
4942 TP c = x[j];
4943 TP c2;
4944 TP coefd[3];
4945 if (c < -1) {
4946 l_border_tmp = -c - 1;
4947 } else if (c > 1) {
4948 l_border_tmp = c - 1;
4949 }
4950 c2 = c * c;
4951 coefd[0] = 0.5 * (c2 - c);
4952 coefd[1] = 0.5 * (c2 + c);
4953 coefd[2] = 1 - c2;
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;
4957 }
4958 if (l_border_tmp < l_border && lr2_tmp < lr2) {
4959 lr2 = lr2_tmp;
4960 l_border = l_border_tmp;
4961 r = c;
4962 }
4963 }
4964 }
4965 TP r2 = r * r;
4966 va = 0.5 * (r2 - r);
4967 vb = 0.5 * (r2 + r);
4968 vc = 1 - r2;
4969 }
4970
4971 // compute the tangent in the quadratic edge defined by a, b and c
4972 // a---------c-----b
4973 // pos = -1 -> a ; pos = 0 -> c ; pos = 1 -> b
4974 void get_tangent(
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) {
4977 tangent[i] =
4978 (a_coor[i] + b_coor[i] - 2 * c_coor[i]) * pos + (b_coor[i] - a_coor[i]) / 2;
4979 }
4980 normalise_3(tangent);
4981 }
4982};
4983
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);
4990} // namespace b2000
4991
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);
5000
5001 const int nb_nodes_of = EDGE::nb_nodes_of;
5002 const int* nodes_of_id = EDGE::nodes_of_id;
5003
5004 // The normal director at the initial configuration for each nodes.
5005 TP D[nb_nodes][3] = {};
5006
5007 // The local referential for the nodes that have 5 dofs.
5008 TP Dr[nb_nodes_of][2][3] = {};
5009
5010 // The coordinate at each nodes
5011 TP R[nb_nodes][6] = {};
5012
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]);
5016 } else {
5017 dof5_node_type::get_coor_s(R[k], nodes[k]);
5018 }
5019 }
5020
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]]) {
5024 number_of_dof += 6;
5025 } else {
5026 number_of_dof += 5;
5027 }
5028 }
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) { // inversed normal w.r.t. element
5034 scale_3(R[k] + 3, -1.);
5035 }
5036 std::copy(R[k] + 3, R[k] + 6, D[k]);
5037 }
5038 }
5039
5040 // update D and R
5041 for (int k = 0; k != nb_nodes_of; ++k) {
5042 const int kk = nodes_of_id[k];
5043 if (kk != k) {
5044 std::copy(D[kk], D[kk + 1], D[k]);
5045 std::copy(R[kk], R[kk + 1], R[k]);
5046 }
5047 }
5048
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};
5052 outer_product_3(e2, D[kk], Dr[kk][0]);
5053 if (normalise_3(Dr[kk][0]) < 0.01) {
5054 TP e2[3] = {1, 0, 0};
5055 outer_product_3(e2, D[kk], Dr[kk][0]);
5056 normalise_3(Dr[kk][0]);
5057 }
5058 outer_product_3(D[kk], Dr[kk][0], Dr[kk][1]);
5059 normalise_3(Dr[kk][1]);
5060 }
5061 }
5062
5063 // The 3 translations at each nodes
5064 TP u_dof[nb_nodes_of][3];
5065
5066 // The 3 rotations at each nodes
5067 TP w_dof[nb_nodes_of][3];
5068
5069 // convert the 5 dof to 6 dof
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);
5073 } else {
5074 size_t o = 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]);
5080 o += 6;
5081 } else {
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];
5086 o += 5;
5087 }
5088 }
5089 }
5090
5091 if (FREEZ) {
5092 // the two vectors orthogonal to D
5093 TP Dn[2][nb_nodes_of][3];
5094 if (nb_nodes_of == 2) {
5095 sub_3(R[0], R[1], Dn[0][0]);
5096 normalise_3(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]);
5102 } else {
5104 }
5105
5106 // The normal to the director in the deformed configuration at each nodes.
5107 TP d[2][nb_nodes_of][3];
5108
5109 // The derivative of the normal to the director in the deformed
5110 // configuration relatively to the rotation dof for each nodes.
5111 TP T[2][nb_nodes_of][3][3];
5112
5113 // Computation of the deformed director and this derivative at
5114 // each nodes. This is done by finite rotation.
5115 for (int kk = 0; kk != nb_nodes_of; ++kk) {
5116 outer_product_3(Dn[0][kk], D[kk], Dn[1][kk]);
5117 normalise_3(Dn[1][kk]);
5118
5119 TP norm_w;
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);
5123 TP O2[6];
5124 {
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;
5134 }
5135 TP s;
5136 TP c;
5137 TP cc;
5138 if (ww < 0.25) {
5139 TP w4 = ww * ww;
5140 TP w6 = w4 * ww;
5141 TP w8 = w6 * ww;
5142 TP w10 = w8 * ww;
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;
5146 } else {
5147 s = std::sin(norm_w) / norm_w;
5148 c = std::sin(norm_w * 0.5) / norm_w;
5149 c = 2 * c * c;
5150 cc = (norm_w - std::sin(norm_w)) / (ww * norm_w);
5151 }
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];
5161 }
5162
5163 if (trans_d_constraint_d_dof.is_null()) { continue; }
5164 s = c;
5165 c = cc;
5166
5167 TP HT3k[3][3];
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];
5178 } else {
5179 {
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];
5187 }
5188 {
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];
5196 }
5197 }
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];
5210 }
5211 }
5212 }
5213
5214 TP c_coor[3] = {};
5215 volume_node_type::get_coor_s(c_coor, nodes[nb_nodes]);
5216 TP v[nb_nodes_of];
5217
5218 if (nb_nodes_of == 2) {
5219 v[1] = 0;
5220 TP l = 0;
5221 for (int i = 0; i != 3; ++i) {
5222 const TP b_a = R[1][i] - R[0][i];
5223 l += b_a * b_a;
5224 v[1] += b_a * (c_coor[i] - R[0][i]);
5225 }
5226 v[1] /= l;
5227 v[0] = 1 - v[1];
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]);
5230 } else {
5232 }
5233
5234 TP cu[3];
5235 TP cd[2][3];
5236 {
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];
5246 }
5247 }
5248 }
5249
5250 if (!constraint.is_null()) {
5251 constraint.resize(2);
5252 if (dof.size2() == 0) {
5253 constraint.set_zero();
5254 } else {
5255 for (int k = 0; k != 2; ++k) {
5256 constraint[k] = 0;
5257 for (int i = 0; i != 3; ++i) { constraint[k] += cu[i] * cd[k][i]; }
5258 }
5259 }
5260 }
5261
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));
5265 size_t* colind;
5266 size_t* rowind;
5267 TP* value;
5268 trans_d_constraint_d_dof.get_values(colind, rowind, value);
5269 size_t* rowind_begin = rowind;
5270 *colind++ = 0;
5271 for (int i = 0; i != 2; ++i) {
5272 int pos = 0;
5273 for (int k = 0; k != nb_nodes_of; ++k) {
5274 TP dd[3];
5275 for (int j = 0; j != 3; ++j, ++pos) {
5276 *rowind++ = pos;
5277 *value++ = -v[k] * d[i][k][j];
5278 dd[j] = cu[j] * v[k];
5279 }
5280 int s = nodes_type_6_dof[nodes_of_id[k]] ? 3 : 2;
5281 for (int j = 0; j != s; ++j, ++pos) {
5282 *rowind++ = pos;
5283 TP tmp = 0;
5284 for (int jj = 0; jj != 3; ++jj) { tmp += dd[jj] * T[i][k][j][jj]; }
5285 *value++ = tmp;
5286 }
5287 }
5288 for (int j = 0; j != 3; ++j, ++pos) {
5289 *rowind++ = pos;
5290 *value++ = cd[i][j];
5291 }
5292 *colind++ = rowind - rowind_begin;
5293 }
5294 }
5295
5296 } else {
5297 // The director in the deformed configuration at each nodes.
5298 TP d[nb_nodes_of][3];
5299
5300 // The derivative of the director in the deformed
5301 // configuration relatively to the rotation dof for each nodes.
5302 TP T[nb_nodes_of][3][3];
5303
5304 // Computation of the deformed director and this derivative at
5305 // each nodes. This is done by finite rotation.
5306 for (int kk = 0; kk != nb_nodes_of; ++kk) {
5307 TP norm_w;
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);
5311 TP O2[6];
5312 {
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;
5322 }
5323 TP s;
5324 TP c;
5325 TP cc;
5326 if (ww < 0.25) {
5327 TP w4 = ww * ww;
5328 TP w6 = w4 * ww;
5329 TP w8 = w6 * ww;
5330 TP w10 = w8 * ww;
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;
5334 } else {
5335 s = std::sin(norm_w) / norm_w;
5336 c = std::sin(norm_w * 0.5) / norm_w;
5337 c = 2 * c * c;
5338 cc = (norm_w - std::sin(norm_w)) / (ww * norm_w);
5339 }
5340 {
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];
5349 }
5350
5351 if (trans_d_constraint_d_dof.is_null()) { continue; }
5352 s = c;
5353 c = cc;
5354
5355 TP HT3k[3][3];
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];
5366 } else {
5367 {
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];
5375 }
5376 {
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];
5384 }
5385 }
5386 {
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];
5398 }
5399 }
5400 }
5401
5402 TP c_coor[3] = {};
5403 volume_node_type::get_coor_s(c_coor, nodes[nb_nodes]);
5404 TP v[nb_nodes_of];
5405
5406 if (nb_nodes_of == 2) {
5407 v[1] = 0;
5408 TP l = 0;
5409 for (int i = 0; i != 3; ++i) {
5410 const TP b_a = R[1][i] - R[0][i];
5411 l += b_a * b_a;
5412 v[1] += b_a * (c_coor[i] - R[0][i]);
5413 }
5414 v[1] /= l;
5415 v[0] = 1 - v[1];
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]);
5418 } else {
5420 }
5421
5422 TP vz = 0;
5423 for (int i = 0; i != 3; ++i) {
5424 TP p = c_coor[i];
5425 TP dir = 0;
5426 for (int j = 0; j != nb_nodes_of; ++j) {
5427 p -= R[j][i] * v[j];
5428 dir += D[j][i] * v[j];
5429 }
5430 vz += p * dir;
5431 }
5432
5433 if (!constraint.is_null()) {
5434 constraint.resize(3);
5435 if (dof.size2() == 0) {
5436 constraint.set_zero();
5437 } else {
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];
5443 }
5444 constraint[i] = ctmp;
5445 }
5446 }
5447 }
5448
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);
5453 size_t* colind;
5454 size_t* rowind;
5455 TP* value;
5456 trans_d_constraint_d_dof.get_values(colind, rowind, value);
5457 size_t* rowind_begin = rowind;
5458 *colind++ = 0;
5459 for (int i = 0; i != 3; ++i) {
5460 int pos = 0;
5461 for (int k = 0; k != nb_nodes_of; ++k) {
5462 *rowind++ = i + pos;
5463 *value++ = v[k];
5464 pos += 3;
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];
5469 }
5470 pos += s;
5471 }
5472 *rowind++ = pos + i;
5473 *value++ = -1;
5474 *colind++ = rowind - rowind_begin;
5475 }
5476 }
5477 }
5478}
5479
5480#ifdef TT_INTEGRATION_AT_GAUSS_POINT
5481#undef TT_INTEGRATION_AT_GAUSS_POINT
5482#endif
5483
5484#endif // B2ELEMENT_SHELL_MITC_H_
#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