b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2reissner_beam_property.H
1//------------------------------------------------------------------------
2// b2reissner_beam_property.H --
3//
4//
5// written by Mathias Doreille, Silvio Merazzi <merazzi@smr.ch>
6// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
7//
8// (c) 2011-2023 SMR Engineering & Development SA
9// 2502 Bienne, Switzerland
10//
11// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
12// Linder Höhe, 51147 Köln
13//
14// All Rights Reserved. Proprietary source code. The contents of
15// this file may not be disclosed to third parties, copied or
16// duplicated in any form, in whole or in part, without the prior
17// written permission of SMR.
18//------------------------------------------------------------------------
19
20// #undef DEBUG
21
22#ifndef B2REISSNER_BEAM_PROPERTY_H_
23#define B2REISSNER_BEAM_PROPERTY_H_
24
25#include <bitset>
26#include <cmath>
27
28#include "b2ppconfig.h"
29#include "b2solid_mechanics.H"
30#include "model/b2element.H"
31#include "model/b2model.H"
32#include "model/b2node.H"
33#include "model/b2solution.H"
34#include "utils/b2beam_cross_sections.H"
35#include "utils/b2exception.H"
36
37namespace b2000 {
38
39class ReissnerBeamPropertyBase : virtual public SolidMechanics {
40public:
41 enum Components {
42 only_membrane_and_shear_components = 1,
43 all_but_membrane_and_shear_components = 2,
44 all_components = 3
45 };
46
47 LinearType linear(const int layer_id = -1) const { return nonlinear; }
48
49 bool isotropic(const int layer_id = -1) const { return false; }
50
51 virtual bool path_dependent(const int layer_id = -1) const { return false; }
52};
53
54template <typename T>
55class ReissnerBeamProperty : public ReissnerBeamPropertyBase {
56 //: virtual public SolidMechanics {
57public:
58 // enum Components {
59 // only_membrane_and_shear_components = 1,
60 // all_but_membrane_and_shear_components = 2,
61 // all_components = 3
62 // };
63
64 // LinearType
65 // linear(const int layer_id = -1) const {
66 // return nonlinear;
67 // }
68
69 // bool
70 // isotropic(const int layer_id = -1) const {
71 // return false;
72 // }
73
74 // virtual bool
75 // path_dependent(const int layer_id = -1) const {
76 // return false;
77 // }
78
102 virtual void get_value(
103 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
104 const double time, const double delta_time, GradientContainer* gradient_container,
105 SolverHints* solver_hints, const Element& element, const double element_coordinate,
106 const b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
107 const T reference_to_undeformed_cross_section_rotator[3][3],
108 const T undeformed_to_deformed_cross_section_rotator[3][3],
109 const T strain_and_curvature[6], const T acceleration[6],
110 const ReissnerBeamPropertyBase::Components components, T force_and_moment[6],
111 T d_force_and_moment_d_strain_and_curvature[6][6], T d_force_and_moment_d_time[6],
112 T inertia[6], T d_inertia_d_acceleration[3][3][3]) = 0;
113
114 virtual bool get_nodes_eccentricity_and_axe(
115 int nb_nodes, const T node_coor[][3], const double node_icoor[], const T node_e3[][3],
116 T e1[3], T nodes_ecc[][3]) {
117 UnimplementedError() << typeid(*this) << " get_nodes_eccentricity_and_axe()" << THROW;
118 return false;
119 }
120
121 virtual const T* get_C() const {
122 UnimplementedError() << typeid(*this) << " get_C()" << THROW;
123 return 0;
124 }
125
126 virtual void set_non_structural_mass(const T& mass_, const T inertia_point_[2]) {}
127
128 virtual void set_pre_stress(const T& force, const T coor[2]) {}
129
130 virtual void set_pre_stress(const T force[3]) {}
131
137 virtual const T* get_cross_section_face_segment(const int face_id) const {
138 UnimplementedError() << " " << typeid(*this) << " get_cross_section_face_segment() "
139 << THROW;
140 return 0;
141 }
142
143 virtual bool get_nodes_neutral_point(
144 int nb_nodes, const T node_coor[][3], const double node_icoor[], const T node_e3[][3],
145 T neutral_point[][2]) {
146 UnimplementedError() << " " << typeid(*this) << " get_nodes_neutral_point() " << THROW;
147 return false;
148 }
149
150 virtual void field_section_integration(
151 const double element_coordinate, const T value[6], const bool multiply_with_density,
152 T force[6]) {
153 UnimplementedError() << " " << typeid(*this) << " field_section_integration() " << THROW;
154 }
155
156 virtual std::bitset<6> get_release_end_node_dof(const int node_ind) const {
157 Exception() << THROW;
158 return std::bitset<6>();
159 }
160
161 T get_prop_required(const RTable rtable, const std::string& s) const {
162 if (rtable.has_key(s)) {
163 return rtable.get<T>(s);
164 } else {
165 Exception() << "Required beam section property " << s << " not found" << THROW;
166 return T(0);
167 }
168 }
169};
170
171template <typename T>
172class ReissnerBeamConstantCrossSectionElasticProperty : public ReissnerBeamProperty<T> {
173public:
174 ReissnerBeamConstantCrossSectionElasticProperty() {
175 std::fill_n(shear_center, 2, T(0));
176 std::fill_n(neutral_axis_point, 2, T(0));
177 std::fill_n(inertia_point, 2, T(0));
178 mass = 0;
179 ns_mass = 0;
180 std::fill_n(ns_inertia_point, 2, T(0));
181 std::fill_n(pre_stress, 3, T(0));
182 area = 0;
183 }
184
185 // void print_array(std::string t, T a[], int n) {
186 // std::cerr << "ReissnerBeam Array " << t;
187 // for (int i = 0; i < n; ++i) std::cerr << " " << a[i];
188 // std::cerr << std::endl;
189 // }
190
191 bool get_nodes_neutral_point(
192 int nb_nodes, const T node_coor[][3], const double node_icoor[], const T node_e3[][3],
193 T neutral_point[][2]) {
194 for (int k = 0; k != nb_nodes; ++k) {
195 neutral_point[k][0] = neutral_axis_point[0];
196 neutral_point[k][1] = neutral_axis_point[1];
197 }
198 return true;
199 }
200
201 const T* get_C() const { return &C[0]; }
202
203 void get_value(
204 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
205 const double time, const double delta_time, GradientContainer* gradient_container,
206 SolverHints* solver_hints, const Element& element, const double element_coordinate,
207 const b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
208 const T reference_to_undeformed_cross_section_rotator[3][3],
209 const T undeformed_to_deformed_cross_section_rotator[3][3],
210 const T strain_and_curvature[6], const T acceleration[6],
211 const ReissnerBeamPropertyBase::Components components, T force_and_moment[6],
212 T d_force_and_moment_d_strain_and_curvature[6][6], T d_force_and_moment_d_time[6],
213 T inertia[6], T d_inertia_d_acceleration[3][3][3]) {
214 const bool neutral_axis_point_zero =
215 true; // neutral_axis_point[0] == 0 && neutral_axis_point[1] == 0;
216
217 if (gradient_container != nullptr && reference_to_undeformed_cross_section_rotator != 0) {
218 static const GradientContainer::FieldDescription descr_mbase = {
219 "MBASE_IP",
220 "e1x.e1y.e1z.e2x.e2y.e2z.e3x.e3y.e3z",
221 "Material base vectors in the undeformed configuration",
222 3,
223 9,
224 1,
225 3,
226 false,
227 typeid(T)};
228 T mbase[3][3];
229 std::copy(
230 reference_to_undeformed_cross_section_rotator[2],
231 reference_to_undeformed_cross_section_rotator[2] + 3, mbase[0]);
232 std::copy(
233 reference_to_undeformed_cross_section_rotator[0],
234 reference_to_undeformed_cross_section_rotator[0] + 3, mbase[1]);
235 std::copy(
236 reference_to_undeformed_cross_section_rotator[1],
237 reference_to_undeformed_cross_section_rotator[1] + 3, mbase[2]);
238 if (gradient_container->compute_field_value(descr_mbase.name)) {
239 gradient_container->set_field_value(descr_mbase, &mbase[0][0]);
240 }
241 }
242
243 if (force_and_moment) {
244 if (strain_and_curvature) {
245 if (components != this->only_membrane_and_shear_components) {
246 for (int i = 3; i != 6; ++i) {
247 force_and_moment[i] = strain_and_curvature[i] * C[i];
248 }
249 force_and_moment[3] += strain_and_curvature[4] * C[7];
250 force_and_moment[4] += strain_and_curvature[3] * C[7];
251 } else {
252 std::fill(force_and_moment + 3, force_and_moment + 6, T(0));
253 }
254 if (components == this->only_membrane_and_shear_components
255 || components == this->all_components) {
256 if (neutral_axis_point_zero) {
257 for (int i = 0; i != 3; ++i) {
258 force_and_moment[i] = C[i] * strain_and_curvature[i];
259 }
260 force_and_moment[0] += strain_and_curvature[1] * C[6];
261 force_and_moment[1] += strain_and_curvature[0] * C[6];
262 } else {
263 force_and_moment[0] =
264 C[0] * strain_and_curvature[0] + strain_and_curvature[1] * C[6];
265 force_and_moment[1] =
266 C[1] * strain_and_curvature[1] + strain_and_curvature[0] * C[6];
267 force_and_moment[2] = C[2]
268 * (strain_and_curvature[2]
269 + neutral_axis_point[1] * strain_and_curvature[3]
270 - neutral_axis_point[0] * strain_and_curvature[4]);
271 force_and_moment[3] += force_and_moment[2] * neutral_axis_point[1];
272 force_and_moment[4] += -force_and_moment[2] * neutral_axis_point[0];
273 }
274 } else {
275 std::fill(force_and_moment, force_and_moment + 3, T(0));
276 }
277 } else {
278 std::fill(force_and_moment, force_and_moment + 6, T(0));
279 }
280 if (components != this->only_membrane_and_shear_components) {
281 for (int i = 0; i != 3; ++i) { force_and_moment[i] += pre_stress[i]; }
282 }
283 }
284
285 if (d_force_and_moment_d_strain_and_curvature) {
286 if (components == this->only_membrane_and_shear_components
287 || components == this->all_components) {
288 std::fill(
289 d_force_and_moment_d_strain_and_curvature[0],
290 d_force_and_moment_d_strain_and_curvature[6], T(0));
291 for (int i = 0; i != 3; ++i) {
292 d_force_and_moment_d_strain_and_curvature[i][i] = C[i];
293 }
294 d_force_and_moment_d_strain_and_curvature[1][0] =
295 d_force_and_moment_d_strain_and_curvature[0][1] = C[6];
296 if (!neutral_axis_point_zero) {
297 d_force_and_moment_d_strain_and_curvature[3][2] =
298 d_force_and_moment_d_strain_and_curvature[2][3] =
299 C[2] * neutral_axis_point[1];
300 d_force_and_moment_d_strain_and_curvature[4][2] =
301 d_force_and_moment_d_strain_and_curvature[2][4] =
302 -C[2] * neutral_axis_point[0];
303 d_force_and_moment_d_strain_and_curvature[3][3] +=
304 C[2] * neutral_axis_point[1] * neutral_axis_point[1];
305 d_force_and_moment_d_strain_and_curvature[5][4] =
306 d_force_and_moment_d_strain_and_curvature[4][5] =
307 -C[2] * neutral_axis_point[0] * neutral_axis_point[1];
308 d_force_and_moment_d_strain_and_curvature[4][4] +=
309 C[2] * neutral_axis_point[0] * neutral_axis_point[0];
310 }
311 } else {
312 std::fill(
313 d_force_and_moment_d_strain_and_curvature[0],
314 d_force_and_moment_d_strain_and_curvature[6], T(0));
315 }
316 if (components != this->only_membrane_and_shear_components) {
317 for (int i = 3; i != 6; ++i) {
318 d_force_and_moment_d_strain_and_curvature[i][i] += C[i];
319 }
320 d_force_and_moment_d_strain_and_curvature[3][4] += C[7];
321 d_force_and_moment_d_strain_and_curvature[4][3] += C[7];
322 }
323 }
324
325 const bool inertia_point_zero = inertia_point[0] == 0 && inertia_point[1] == 0;
326 const bool ns_inertia_point_zero = ns_inertia_point[0] == 0 && ns_inertia_point[1] == 0;
327 if (inertia) {
328 if (acceleration) {
329 if (inertia_point_zero) {
330 for (int i = 0; i != 3; ++i) { inertia[i] = mass * acceleration[i]; }
331 for (int i = 0; i != 2; ++i) {
332 inertia[i + 3] = inertia_moment[i] * acceleration[i + 3];
333 }
334 inertia[5] = (inertia_moment[0] + inertia_moment[1]) * acceleration[5];
335 inertia[3] += inertia_moment[2] * acceleration[4];
336 inertia[4] += inertia_moment[2] * acceleration[3];
337 } else {
338 inertia[0] = mass * (acceleration[0] - inertia_point[1] * acceleration[5]);
339 inertia[1] = mass * (acceleration[1] + inertia_point[0] * acceleration[5]);
340 inertia[2] = mass
341 * (acceleration[2] + inertia_point[1] * acceleration[3]
342 - inertia_point[0] * acceleration[4]);
343 inertia[3] = inertia_moment[0] * acceleration[3]
344 + inertia_moment[2] * acceleration[4]
345 + inertia[2] * inertia_point[1];
346 inertia[4] = inertia_moment[1] * acceleration[4]
347 + inertia_moment[2] * acceleration[3]
348 - inertia[2] * inertia_point[0];
349 inertia[5] = (inertia_moment[0] + inertia_moment[1]) * acceleration[5]
350 - inertia[0] * inertia_point[1] + inertia[1] * inertia_point[1];
351 }
352 if (ns_inertia_point_zero) {
353 for (int i = 0; i != 3; ++i) { inertia[i] += ns_mass * acceleration[i]; }
354 } else {
355 const T i0 =
356 ns_mass * (acceleration[0] - ns_inertia_point[1] * acceleration[5]);
357 inertia[0] += i0;
358 const T i1 =
359 ns_mass * (acceleration[1] + ns_inertia_point[0] * acceleration[5]);
360 inertia[1] += i1;
361 const T i2 = ns_mass
362 * (acceleration[2] + ns_inertia_point[1] * acceleration[3]
363 - ns_inertia_point[0] * acceleration[4]);
364 inertia[2] += i2;
365 inertia[3] += i2 * ns_inertia_point[1];
366 inertia[4] += -i2 * ns_inertia_point[0];
367 inertia[5] += -i0 * ns_inertia_point[1] + i1 * ns_inertia_point[1];
368 }
369 } else {
370 std::fill(inertia, inertia + 6, T(0));
371 }
372 }
373
374 if (d_inertia_d_acceleration) {
375 std::fill(d_inertia_d_acceleration[0][0], d_inertia_d_acceleration[3][0], T(0));
376 for (int i = 0; i != 3; ++i) { d_inertia_d_acceleration[0][i][i] = mass + ns_mass; }
377 for (int i = 0; i != 2; ++i) { d_inertia_d_acceleration[2][i][i] = inertia_moment[i]; }
378 d_inertia_d_acceleration[2][2][2] = inertia_moment[0] + inertia_moment[1];
379 d_inertia_d_acceleration[2][1][0] = inertia_moment[2];
380 if (!inertia_point_zero) {
381 d_inertia_d_acceleration[1][0][2] =
382 mass * inertia_point[1] + ns_mass * ns_inertia_point[1];
383 d_inertia_d_acceleration[1][1][2] =
384 -mass * inertia_point[0] - ns_mass * ns_inertia_point[0];
385 d_inertia_d_acceleration[1][2][0] =
386 -mass * inertia_point[1] - ns_mass * ns_inertia_point[1];
387 d_inertia_d_acceleration[1][2][1] =
388 mass * inertia_point[0] + ns_mass * ns_inertia_point[0];
389 d_inertia_d_acceleration[2][0][0] +=
390 mass * inertia_point[1] * inertia_point[1]
391 + ns_mass * ns_inertia_point[1] * ns_inertia_point[1];
392 d_inertia_d_acceleration[2][1][0] +=
393 -mass * inertia_point[0] * inertia_point[1]
394 - -ns_mass * ns_inertia_point[0] * ns_inertia_point[1];
395 d_inertia_d_acceleration[2][1][1] +=
396 mass * inertia_point[0] * inertia_point[0]
397 + ns_mass * ns_inertia_point[0] * ns_inertia_point[0];
398 d_inertia_d_acceleration[2][2][2] +=
399 mass
400 * (inertia_point[0] * inertia_point[0]
401 + inertia_point[1] * inertia_point[1])
402 + ns_mass
403 * (ns_inertia_point[0] * ns_inertia_point[0]
404 + ns_inertia_point[1] * ns_inertia_point[1]);
405 }
406 d_inertia_d_acceleration[2][0][1] = d_inertia_d_acceleration[2][1][0];
407 }
408
409 if (d_force_and_moment_d_time) {
410 std::fill(d_force_and_moment_d_time, d_force_and_moment_d_time + 6, T(0));
411 }
412
413 if (gradient_container
414 && (components == this->only_membrane_and_shear_components
415 || components == this->all_components)) {
416 static const GradientContainer::FieldDescription force_beam_descr = {
417 "FORCE_SECTION_BEAM",
418 "Fx.Fy.Fz",
419 "Beam section forces",
420 3,
421 3,
422 1,
423 3,
424 false,
425 typeid(T)};
426 if (gradient_container->compute_field_value(force_beam_descr.name)) {
427 const T force[3] = {
428 C[2] * strain_and_curvature[2],
429 C[0] * strain_and_curvature[0] + strain_and_curvature[1] * C[6],
430 C[1] * strain_and_curvature[1] + strain_and_curvature[0] * C[6]};
431 gradient_container->set_field_value(force_beam_descr, force);
432 }
433 static const GradientContainer::FieldDescription moment_beam_descr = {
434 "MOMENT_SECTION_BEAM",
435 "Mx.My.Mz",
436 "Beam section moments",
437 3,
438 3,
439 1,
440 3,
441 false,
442 typeid(T)};
443 if (gradient_container->compute_field_value(moment_beam_descr.name)) {
444 const T moment[3] = {
445 C[5] * strain_and_curvature[5],
446 C[3] * strain_and_curvature[3] + strain_and_curvature[4] * C[7],
447 C[4] * strain_and_curvature[4] + strain_and_curvature[3] * C[7]};
448 gradient_container->set_field_value(moment_beam_descr, moment);
449 }
450
451 set_gradient(
452 reference_to_undeformed_cross_section_rotator,
453 undeformed_to_deformed_cross_section_rotator, strain_and_curvature,
454 gradient_container);
455 }
456 }
457
458 void field_section_integration(
459 const double element_coordinate, const T value[6], const bool multiply_with_density,
460 T force[6]) {
461 if (multiply_with_density) {
462 // Mass.
463 force[0] = mass * (value[0] - inertia_point[1] * value[5]);
464 force[1] = mass * (value[1] + inertia_point[0] * value[5]);
465 force[2] =
466 mass * (value[2] + inertia_point[1] * value[3] - inertia_point[0] * value[4]);
467 force[3] = inertia_moment[0] * value[3] + inertia_moment[2] * value[4]
468 + force[2] * inertia_point[1];
469 force[4] = inertia_moment[1] * value[4] + inertia_moment[2] * value[3]
470 - value[2] * inertia_point[0];
471 force[5] = (inertia_moment[0] + inertia_moment[1]) * value[5]
472 - force[0] * inertia_point[1] + force[1] * inertia_point[1];
473
474 // Contribution of nonstructural mass.
475 const T i0 = ns_mass * (value[0] - ns_inertia_point[1] * value[5]);
476 force[0] += i0;
477 const T i1 = ns_mass * (value[1] + ns_inertia_point[0] * value[5]);
478 force[1] += i1;
479 const T i2 =
480 ns_mass
481 * (value[2] + ns_inertia_point[1] * value[3] - ns_inertia_point[0] * value[4]);
482 force[2] += i2;
483 force[3] += i2 * ns_inertia_point[1];
484 force[4] += -i2 * ns_inertia_point[0];
485 force[5] += -i0 * ns_inertia_point[1] + i1 * ns_inertia_point[1];
486 } else {
487 force[0] = area * value[0];
488 force[1] = area * value[1];
489 force[2] = area * value[2];
490 }
491 }
492
493 virtual void set_gradient(
494 const T reference_to_undeformed_cross_section_rotator[3][3],
495 const T undeformed_to_deformed_cross_section_rotator[3][3],
496 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
497
498 void set_non_structural_mass(const T& mass_, const T inertia_point_[2]) {
499 ns_mass = mass_;
500 std::copy(inertia_point_, inertia_point_ + 2, ns_inertia_point);
501 }
502
503 void set_pre_stress(const T& force, const T coor[2]) {
504 pre_stress[2] = force;
505 if (coor) {
506 pre_stress[0] = (coor[1] - neutral_axis_point[1]) * force;
507 pre_stress[1] = -(coor[0] - neutral_axis_point[0]) * force;
508 }
509 }
510
511 void set_pre_stress(const T force[3]) {
512 pre_stress[2] = force[0];
513 pre_stress[0] = force[1];
514 pre_stress[1] = force[2];
515 }
516
517public:
518 T shear_center[2];
519 T neutral_axis_point[2];
520 T C[8];
521 T inertia_point[2];
522 T mass;
523 T inertia_moment[3];
524 T ns_inertia_point[2];
525 T ns_mass;
526 T pre_stress[3];
527 T area;
528};
529
530template <typename T>
531class ReissnerBeamVariableCrossSectionElasticProperty : public ReissnerBeamProperty<T> {
532public:
533 bool get_nodes_neutral_point(
534 int nb_nodes, const T node_coor[][3], const double node_icoor[], const T node_e3[][3],
535 T neutral_point[][2]) {
536 for (int k = 0; k != nb_nodes; ++k) {
537 const size_t i1 =
538 std::lower_bound(property_icoor.begin(), property_icoor.end(), node_icoor[k])
539 - property_icoor.begin();
540 const size_t i2 = i1 + 1;
541 const T s1 =
542 (node_icoor[k] - property_icoor[i2]) / (property_icoor[i1] - property_icoor[i2]);
543 const T s2 = 1 - s1;
544 for (int i = 0; i != 2; ++i) {
545 neutral_point[k][i] =
546 s1 * property[i1].shear_center[i] + s2 * property[i2].shear_center[i];
547 }
548 }
549 return true;
550 }
551
552 void get_value(
553 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
554 const double time, const double delta_time, const Element& element,
555 const double element_coordinate,
556 const b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
557 const T reference_to_undeformed_cross_section_rotator[3][3],
558 const T undeformed_to_deformed_cross_section_rotator[3][3],
559 const T strain_and_curvature[6], const T acceleration[6],
560 const ReissnerBeamPropertyBase::Components components, T force_and_moment[6],
561 T d_force_and_moment_d_strain_and_curvature[6][6], T d_force_and_moment_d_time[6],
562 T inertia[6], T d_inertia_d_acceleration[3][3][3], GradientContainer* gradient_container,
563 SolverHints* solver_hints) {
564 const size_t i1 =
565 std::lower_bound(property_icoor.begin(), property_icoor.end(), element_coordinate)
566 - property_icoor.begin();
567 const size_t i2 = i1 + 1;
568 const T s1 =
569 (element_coordinate - property_icoor[i2]) / (property_icoor[i1] - property_icoor[i2]);
570 const T s2 = 1 - s1;
571
572 ReissnerBeamConstantCrossSectionElasticProperty<T> cs;
573 for (int i = 0; i != 2; ++i) {
574 cs.neutral_axis_point[i] =
575 s1 * property[i1].neutral_axis_point[i] + s2 * property[i2].neutral_axis_point[i];
576 cs.inertia_point[i] =
577 s1 * property[i1].inertia_point[i] + s2 * property[i2].inertia_point[i];
578 cs.ns_inertia_point[i] =
579 s1 * property[i1].ns_inertia_point[i] + s2 * property[i2].ns_inertia_point[i];
580 }
581 for (int i = 0; i != 7; ++i) { cs.C[i] = s1 * property[i1].C[i] + s2 * property[i2].C[i]; }
582 cs.mass = s1 * property[i1].mass + s2 * property[i2].mass;
583 cs.ns_mass = s1 * property[i1].ns_mass + s2 * property[i2].ns_mass;
584 for (int i = 0; i != 3; ++i) {
585 cs.inertia_moment[i] =
586 s1 * property[i1].inertia_moment[i] + s2 * property[i2].inertia_moment[i];
587 }
588 for (int i = 0; i != 3; ++i) {
589 cs.pre_stress[i] = s1 * property[i1].pre_stress[i] + s2 * property[i2].pre_stress[i];
590 }
591
592 cs.get_value(
593 model, linear, equilibrium_solution, time, delta_time, gradient_container,
594 solver_hints, element, element_coordinate, nodes_interpolation,
595 reference_to_undeformed_cross_section_rotator,
596 undeformed_to_deformed_cross_section_rotator, strain_and_curvature, acceleration,
597 components, force_and_moment, d_force_and_moment_d_strain_and_curvature,
598 d_force_and_moment_d_time, inertia, d_inertia_d_acceleration);
599 }
600
601 void field_section_integration(
602 const T element_coordinate, const T value[6], const bool multiply_with_density,
603 T force[6]) {
604 const size_t i1 =
605 std::lower_bound(property_icoor.begin(), property_icoor.end(), element_coordinate)
606 - property_icoor.begin();
607 const size_t i2 = i1 + 1;
608 const T s1 =
609 (element_coordinate - property_icoor[i2]) / (property_icoor[i1] - property_icoor[i2]);
610 const T s2 = 1 - s1;
611
612 ReissnerBeamConstantCrossSectionElasticProperty<T> cs;
613 for (int i = 0; i != 2; ++i) {
614 cs.neutral_axis_point[i] =
615 s1 * property[i1].neutral_axis_point[i] + s2 * property[i2].neutral_axis_point[i];
616 cs.inertia_point[i] =
617 s1 * property[i1].inertia_point[i] + s2 * property[i2].inertia_point[i];
618 cs.ns_inertia_point[i] =
619 s1 * property[i1].ns_inertia_point[i] + s2 * property[i2].ns_inertia_point[i];
620 }
621 for (int i = 0; i != 7; ++i) { cs.C[i] = s1 * property[i1].C[i] + s2 * property[i2].C[i]; }
622 cs.mass = s1 * property[i1].mass + s2 * property[i2].mass;
623 cs.ns_mass = s1 * property[i1].ns_mass + s2 * property[i2].ns_mass;
624 for (int i = 0; i != 3; ++i) {
625 cs.inertia_moment[i] =
626 s1 * property[i1].inertia_moment[i] + s2 * property[i2].inertia_moment[i];
627 }
628 for (int i = 0; i != 3; ++i) {
629 cs.pre_stress[i] = s1 * property[i1].pre_stress[i] + s2 * property[i2].pre_stress[i];
630 }
631
632 cs.field_section_integration(element_coordinate, value, multiply_with_density, force);
633 }
634
635protected:
636 struct CrossSectionProperty {
637 T shear_center[2];
638 T neutral_axis_point[2];
639 T C[7];
640 T inertia_point[2];
641 T mass;
642 T inertia_moment[3];
643 T ns_inertia_point[2];
644 T ns_mass;
645 T pre_stress[3];
646 };
647 std::vector<T> property_icoor;
648 std::vector<CrossSectionProperty> property;
649};
650
651template <typename T>
652class ReissnerBeamElasticConstantProperty
653 : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
654public:
655 ReissnerBeamElasticConstantProperty() {
656 this->mass = this->ns_mass = 0;
657 std::fill_n(this->inertia_point, 2, T(0));
658 std::fill_n(this->inertia_moment, 3, T(0));
659 std::fill_n(this->ns_inertia_point, 2, T(0));
660 std::fill_n(this->pre_stress, 3, T(0));
661 this->area = 0;
662 }
663
664 const T* get_C() const { return &this->C[0]; }
665
666 void print_constants() {
667 std::cerr << "ReissnerBeamElasticConstantPropertyCross section " << std::endl
668 << "inertia_point=" << this->inertia_point[0] << " " << this->inertia_point[1]
669 << std::endl
670 << "shear_center=" << this->shear_center[0] << " " << this->shear_center[1]
671 << std::endl;
672 }
673
674 void set_static_property(
675 const T neutral_axis_point_[2], const T shear_center_[2], T traction, const T bending[3],
676 T torsion, const T shear[3]) {
677 this->shear_center[0] = shear_center_[0];
678 this->shear_center[1] = shear_center_[1];
679 this->neutral_axis_point[0] = neutral_axis_point_[0];
680 this->neutral_axis_point[1] = neutral_axis_point_[1];
681 this->C[0] = shear[0];
682 this->C[1] = shear[1];
683 this->C[2] = traction;
684 this->C[3] = bending[0];
685 this->C[4] = bending[1];
686 this->C[5] = torsion;
687 this->C[6] = shear[2];
688 this->C[7] = bending[2];
689 // #ifdef DEBUG
690 // print_constants();
691 // this->print_array("ReissnerBeamElasticConstantProperty C=", this->C, 8);
692 // this->print_array("neutral_axis_point=", this->neutral_axis_point, 2);
693 // this->print_array("shear_center=", this->shear_center, 2);
694 // #endif
695 }
696
697 void set_dynamic_property(
698 const T mass_, const T mass_inertia_point_[2], const T mass_inertia_moment_[3],
699 T ns_mass_, const T ns_inertia_point_[2]) {
700 this->mass = mass_;
701 this->inertia_point[0] = mass_inertia_point_[0];
702 this->inertia_point[1] = mass_inertia_point_[1];
703 this->inertia_moment[0] = mass_inertia_moment_[0];
704 this->inertia_moment[1] = mass_inertia_moment_[1];
705 this->inertia_moment[2] = mass_inertia_moment_[2];
706 this->ns_mass = ns_mass_;
707 this->ns_inertia_point[0] = ns_inertia_point_[0];
708 this->ns_inertia_point[1] = ns_inertia_point_[1];
709 }
710
711 void set_gradient(
712 const T reference_to_undeformed_cross_section_rotator[3][3],
713 const T undeformed_to_deformed_cross_section_rotator[3][3],
714 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
715};
716
717//
718// Start BEAM SECTION PROPERTIES
719//
720
721template <typename T>
722class ReissnerBeamElasticBarProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
723public:
724 T set_property(
725 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
726 const RTable& frtable) {
727 const T G = E / (2 * (1 + p));
728 // const T K = 10 * (1 + p) / (12 + 11 * p);
729
730 BeamCrossSectionBar<T> section;
731 section.load(rtable);
732 this->C[0] = G * section.A * section.K1;
733 this->C[1] = G * section.A * section.K2;
734 this->C[2] = E * section.A;
735 this->C[3] = E * section.Iyy;
736 this->C[4] = E * section.Izz;
737 this->C[5] = G * section.It;
738 this->C[6] = 0;
739 this->C[7] = 0;
740 this->mass = density * section.A;
741 this->ns_mass = nsmass;
742 this->inertia_moment[0] = density * section.Iyy;
743 this->inertia_moment[1] = density * section.Izz;
744 this->inertia_moment[2] = 0;
745 this->area = section.A;
746
747 return this->area;
748 }
749
750 const T* get_C() const { return &this->C[0]; }
751
752 void set_gradient(
753 const T reference_to_undeformed_cross_section_rotator[3][3],
754 const T undeformed_to_deformed_cross_section_rotator[3][3],
755 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
756};
757
758template <typename T>
759class ReissnerBeamElasticBoxProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
760public:
761 T set_property(
762 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
763 const RTable& frtable) {
764 // const T K = 20 * (1 + p) / (48 + 39 * p);
765 BeamCrossSectionBox<T> section;
766 section.load(rtable);
767 const T G = E / (2 * (1 + p));
768 this->C[0] = G * section.A * section.K1; // K;
769 this->C[1] = G * section.A * section.K2; // C[0];
770 this->C[2] = E * section.A;
771 this->C[3] = E * section.Iyy;
772 this->C[4] = E * section.Izz;
773 this->C[5] = G * section.It;
774 this->C[6] = 0;
775 this->C[7] = 0;
776 this->mass = density * section.A;
777 this->ns_mass = nsmass;
778 this->inertia_moment[0] = density * section.Iyy;
779 this->inertia_moment[1] = density * section.Izz;
780 this->inertia_moment[2] = 0;
781 this->area = section.A;
782 return this->area;
783 }
784
785 const T* get_C() { return &this->C[0]; }
786
787 void set_gradient(
788 const T reference_to_undeformed_cross_section_rotator[3][3],
789 const T undeformed_to_deformed_cross_section_rotator[3][3],
790 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
791};
792
793template <typename T>
794class ReissnerBeamElasticTubeProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
795public:
796 T set_property(
797 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
798 const RTable& frtable) {
799 BeamCrossSectionTube<T> section;
800 section.load(rtable);
801 const T G = E / (2 * (1 + p));
802 this->C[0] = G * section.A * section.K1;
803 this->C[1] = this->C[0];
804 this->C[2] = E * section.A;
805 this->C[3] = E * section.Iyy;
806 this->C[4] = this->C[3];
807 this->C[5] = G * section.It; // Torsion
808 this->C[6] = 0;
809 this->C[7] = 0;
810 this->mass = density * section.A;
811 this->ns_mass = nsmass;
812 this->inertia_moment[0] = density * section.Iyy;
813 this->inertia_moment[1] = this->inertia_moment[0];
814 this->inertia_moment[2] = 0;
815 this->area = section.A;
816 return section.A;
817 }
818
819 const T* get_C() const { return &this->C[0]; }
820
821 void set_gradient(
822 const T reference_to_undeformed_cross_section_rotator[3][3],
823 const T undeformed_to_deformed_cross_section_rotator[3][3],
824 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
825};
826
827template <typename T>
828class ReissnerBeamElasticLProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
829public:
830 T set_property(
831 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
832 const RTable& frtable) {
833 BeamCrossSectionL<T> section;
834 section.load(rtable);
835 this->shear_center[0] = -section.yc;
836 this->shear_center[1] = -section.zc;
837 const T G = E / (2 * (1 + p));
838 this->C[0] = G * section.A * section.K1;
839 this->C[1] = G * section.A * section.K2;
840 this->C[2] = E * section.A;
841 this->C[3] = E * section.Iyy;
842 this->C[4] = E * section.Izz;
843 this->C[5] = G * section.It; // Torsion
844 this->C[6] = 0;
845 this->C[7] = E * section.Iyz;
846 this->mass = density * section.A;
847 this->ns_mass = nsmass;
848 this->inertia_moment[0] = density * section.Iyy;
849 this->inertia_moment[1] = density * section.Izz;
850 this->inertia_moment[2] = 0;
851 this->area = section.A;
852 return this->area;
853 }
854
855 const T* get_C() const { return &this->C[0]; }
856
857 void set_gradient(
858 const T reference_to_undeformed_cross_section_rotator[3][3],
859 const T undeformed_to_deformed_cross_section_rotator[3][3],
860 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
861};
862
865template <typename T>
866class ReissnerBeamElasticIProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
867public:
868 T set_property(
869 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
870 const RTable& frtable) {
871 BeamCrossSectionI<T> section;
872 section.load(rtable);
873
874 // T d = web_height - 2 * flange_thickness;
875 // const T Iyy = (flange_width * std::pow((d + 2 * flange_thickness), 3) - (flange_width -
876 // web_thickness) * std::pow(d, 3)) / 12; const T Izz = (std::pow(flange_width, 3) *
877 // flange_thickness) / 6. + (std::pow(web_thickness, 3) * d) / 12; const T A = 2. *
878 // flange_width * flange_thickness + web_thickness * d; const T beta = 10 * (1 + p) / (12 +
879 // 11 * p); const T A1 = 2 * flange_width * flange_thickness * beta; const T A2 = web_height
880 // * web_thickness * beta;
881 // // Roark page 401 ff, 800 and neglected 1 term
882 // const T D = (std::pow(flange_thickness, 2) + 0.25 * std::pow(web_thickness,
883 // 2))/flange_thickness; const T K1 = 0.333333333333333 * flange_width *
884 // std::pow(flange_thickness, 3); const T K2 = 0.333333333333333 * (d - 2 *
885 // flange_thickness) * std::pow(web_thickness, 3); T t, t1; if (flange_thickness <
886 // web_thickness)
887 // t = flange_thickness;
888 // else
889 // t = web_thickness;
890 // if (flange_thickness > web_thickness)
891 // t1= flange_thickness;
892 // else
893 // t1 = web_thickness;
894
895 // const T alfa = 0.15 * t / t1;
896 // const T K = 2.0 * section.K1 + section.K2 + 2.0 * alfa * std::pow(D, 4);
897 const T G = E / (2 * (1 + p));
898 this->C[0] = G * section.A * section.K1;
899 this->C[1] = G * section.A * section.K2;
900 this->C[2] = E * section.A;
901 this->C[3] = E * section.Iyy;
902 this->C[4] = E * section.Izz;
903 this->C[5] = G * section.It; // Torsion
904 this->C[6] = 0;
905 this->C[7] = 0; // E * section.Izz;
906 this->mass = density * section.A;
907 this->ns_mass = nsmass;
908 this->inertia_moment[0] = density * section.Iyy;
909 this->inertia_moment[1] = density * section.Izz;
910 this->inertia_moment[2] = 0;
911 this->area = section.A;
912 return section.A;
913 }
914
915 const T* get_C() const { return &this->C[0]; }
916};
917
920template <typename T>
921class ReissnerBeamElasticCProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
922public:
923 T set_property(
924 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
925 const RTable& frtable) {
926 BeamCrossSectionC<T> section;
927 section.load(rtable);
928 this->shear_center[0] = section.ys - section.yc;
929 this->shear_center[1] = section.zs - section.zc;
930 this->neutral_axis_point[0] = 0;
931 this->neutral_axis_point[1] = 0;
932 // T Iyy = section.Iyy; //+ pow(section.zs+section.zc,2)*section.A;
933 const T G = E / (2 * (1 + p));
934 this->C[0] = G * section.A * section.K1;
935 this->C[1] = G * section.A * section.K2;
936 this->C[2] = E * section.A;
937 this->C[3] = E * section.Iyy;
938 this->C[4] = E * section.Izz;
939 this->C[5] = G * section.It; // Torsion
940 this->C[6] = 0;
941 this->C[7] = 0; // E * section.Izz;
942 this->mass = density * section.A;
943 this->ns_mass = nsmass;
944 this->inertia_moment[0] = density * section.Iyy;
945 this->inertia_moment[1] = density * section.Izz;
946 this->inertia_moment[2] = 0;
947 this->area = section.A;
948 return section.A;
949 }
950
951 const T* get_C() const { return &this->C[0]; }
952
953 void set_gradient(
954 const T reference_to_undeformed_cross_section_rotator[3][3],
955 const T undeformed_to_deformed_cross_section_rotator[3][3],
956 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
957};
958
959template <typename T>
960class ReissnerBeamElasticRodProperty : public ReissnerBeamConstantCrossSectionElasticProperty<T> {
961public:
962 T set_property(
963 const T& E, const T& p, const T& density, const T& nsmass, const RTable& rtable,
964 const RTable& frtable) {
965 const T G = E / (2 * (1 + p));
966 // const T K = 10 * (1 + p) / (12 + 11 * p);
967
968 BeamCrossSectionRod<T> section;
969 section.load(rtable);
970 this->C[0] = G * section.A * section.K1;
971 this->C[1] = G * section.A * section.K2;
972 this->C[2] = E * section.A;
973 this->C[3] = E * section.Iyy;
974 this->C[4] = E * section.Izz;
975 this->C[5] = G * section.It;
976 this->C[6] = 0;
977 this->C[7] = 0;
978 this->mass = density * section.A;
979 this->ns_mass = nsmass;
980 this->inertia_moment[0] = density * section.Iyy;
981 this->inertia_moment[1] = density * section.Izz;
982 this->inertia_moment[2] = 0;
983 this->area = section.A;
984 return this->area;
985 }
986
987 const T* get_C() const { return &this->C[0]; }
988
989 void set_gradient(
990 const T reference_to_undeformed_cross_section_rotator[3][3],
991 const T undeformed_to_deformed_cross_section_rotator[3][3],
992 const T strain_and_curvature[6], GradientContainer* gradient_container) {}
993};
994
995template <typename T>
996class ReissnerBeamConstantCrossSectionElasticFEMProperty : public ReissnerBeamProperty<T> {
997public:
998 ReissnerBeamConstantCrossSectionElasticFEMProperty(Model* m)
999 : cs_model(m), mass(0), ns_mass(0), compute_gradients(false) {
1000 std::fill_n(neutral_point, 2, T(0));
1001 std::fill_n(inertia_point, 2, T(0));
1002 std::fill_n(inertia_moment, 3, T(0));
1003 std::fill_n(ns_inertia_point, 2, T(0));
1004 std::fill_n(constant_force_and_moment, 6, T(0));
1005 }
1006
1007 ReissnerBeamConstantCrossSectionElasticFEMProperty()
1008 : cs_model(nullptr), mass(0), ns_mass(0), compute_gradients(false) {
1009 std::fill_n(neutral_point, 2, T(0));
1010 std::fill_n(inertia_point, 2, T(0));
1011 std::fill_n(inertia_moment, 3, T(0));
1012 std::fill_n(ns_inertia_point, 2, T(0));
1013 std::fill_n(constant_force_and_moment, 6, T(0));
1014 }
1015
1016 ~ReissnerBeamConstantCrossSectionElasticFEMProperty() { delete cs_model; }
1017
1018 bool get_nodes_neutral_point(
1019 int nb_nodes, const T node_coor[][3], const T node_icoor[], const T node_e3[][3],
1020 T neutral_point_[][2]) {
1021 for (int k = 0; k != nb_nodes; ++k) {
1022 neutral_point_[k][0] = neutral_point[0];
1023 neutral_point_[k][1] = neutral_point[1];
1024 }
1025 return true;
1026 }
1027
1028 void get_value(
1029 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
1030 const double time, const double delta_time, GradientContainer* gradient_container,
1031 SolverHints* solver_hints, const Element& element, const double element_coordinate,
1032 const b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
1033 const T reference_to_undeformed_cross_section_rotator[3][3],
1034 const T undeformed_to_deformed_cross_section_rotator[3][3],
1035 const T strain_and_curvature[6], const T acceleration[6],
1036 const ReissnerBeamPropertyBase::Components components, T force_and_moment[6],
1037 T d_force_and_moment_d_strain_and_curvature[6][6], T d_force_and_moment_d_time[6],
1038 T inertia[6], T d_inertia_d_acceleration[3][3][3]) {
1039 if (force_and_moment) {
1040 if (strain_and_curvature) {
1041 if (components == this->only_membrane_and_shear_components
1042 || components == this->all_components) {
1043 for (int i = 0; i != 3; ++i) {
1044 T tmp(0);
1045 for (int j = 0; j != 3; ++j) {
1046 tmp += lambda(i, j) * strain_and_curvature[j];
1047 }
1048 force_and_moment[i] = tmp;
1049 }
1050 } else {
1051 std::fill(force_and_moment, force_and_moment + 3, T(0));
1052 }
1053 if (components != this->only_membrane_and_shear_components) {
1054 for (int i = 0; i != 3; ++i) {
1055 T tmp(0);
1056 for (int j = 3; j != 6; ++j) {
1057 tmp += lambda(i, j) * strain_and_curvature[j];
1058 }
1059 force_and_moment[i] += tmp;
1060 }
1061 for (int i = 3; i != 6; ++i) {
1062 T tmp(0);
1063 for (int j = 0; j != 6; ++j) {
1064 tmp += lambda(i, j) * strain_and_curvature[j];
1065 }
1066 force_and_moment[i] = tmp;
1067 }
1068 } else {
1069 std::fill(force_and_moment + 3, force_and_moment + 6, T(0));
1070 }
1071 } else {
1072 std::fill(force_and_moment, force_and_moment + 6, T(0));
1073 }
1074 if (components != this->only_membrane_and_shear_components) {
1075 for (int i = 0; i != 6; ++i) {
1076 force_and_moment[i] += constant_force_and_moment[i];
1077 }
1078 }
1079 }
1080
1081 if (d_force_and_moment_d_strain_and_curvature) {
1082 if (components == this->all_components) {
1083 std::copy(
1084 &lambda(0, 0), &lambda(0, 7), d_force_and_moment_d_strain_and_curvature[0]);
1085 } else {
1086 std::fill(
1087 d_force_and_moment_d_strain_and_curvature[0],
1088 d_force_and_moment_d_strain_and_curvature[6], T(0));
1089 if (components == this->only_membrane_and_shear_components) {
1090 for (int j = 0; j != 3; ++j) {
1091 for (int i = 0; i != 3; ++i) {
1092 d_force_and_moment_d_strain_and_curvature[j][i] = lambda(i, j);
1093 }
1094 }
1095 } else if (components != this->only_membrane_and_shear_components) {
1096 for (int j = 0; j != 3; ++j) {
1097 for (int i = 3; i != 6; ++i) {
1098 d_force_and_moment_d_strain_and_curvature[j][i] = lambda(i, j);
1099 }
1100 }
1101 for (int j = 3; j != 6; ++j) {
1102 for (int i = 0; i != 6; ++i) {
1103 d_force_and_moment_d_strain_and_curvature[j][i] = lambda(i, j);
1104 }
1105 }
1106 }
1107 }
1108 }
1109
1110 const bool inertia_point_zero = inertia_point[0] == 0 && inertia_point[1] == 0;
1111 const bool ns_inertia_point_zero = ns_inertia_point[0] == 0 && ns_inertia_point[1] == 0;
1112 if (inertia) {
1113 if (acceleration) {
1114 if (inertia_point_zero) {
1115 for (int i = 0; i != 3; ++i) { inertia[i] = mass * acceleration[i]; }
1116 for (int i = 0; i != 2; ++i) {
1117 inertia[i + 3] = inertia_moment[i] * acceleration[i + 3];
1118 }
1119 inertia[5] = (inertia_moment[0] + inertia_moment[1]) * acceleration[5];
1120 inertia[3] += inertia_moment[2] * acceleration[4];
1121 inertia[4] += inertia_moment[2] * acceleration[3];
1122 } else {
1123 inertia[0] = mass * (acceleration[0] - inertia_point[1] * acceleration[5]);
1124 inertia[1] = mass * (acceleration[1] + inertia_point[0] * acceleration[5]);
1125 inertia[2] = mass
1126 * (acceleration[2] + inertia_point[1] * acceleration[3]
1127 - inertia_point[0] * acceleration[4]);
1128 inertia[3] = inertia_moment[0] * acceleration[3]
1129 + inertia_moment[2] * acceleration[4]
1130 + inertia[2] * inertia_point[1];
1131 inertia[4] = inertia_moment[1] * acceleration[4]
1132 + inertia_moment[2] * acceleration[3]
1133 - inertia[2] * inertia_point[0];
1134 inertia[5] = (inertia_moment[0] + inertia_moment[1]) * acceleration[5]
1135 - inertia[0] * inertia_point[1] + inertia[1] * inertia_point[1];
1136 }
1137 if (ns_inertia_point_zero) {
1138 for (int i = 0; i != 3; ++i) { inertia[i] += ns_mass * acceleration[i]; }
1139 } else {
1140 const T i0 =
1141 ns_mass * (acceleration[0] - ns_inertia_point[1] * acceleration[5]);
1142 inertia[0] += i0;
1143 const T i1 =
1144 ns_mass * (acceleration[1] + ns_inertia_point[0] * acceleration[5]);
1145 inertia[1] += i1;
1146 const T i2 = ns_mass
1147 * (acceleration[2] + ns_inertia_point[1] * acceleration[3]
1148 - ns_inertia_point[0] * acceleration[4]);
1149 inertia[2] += i2;
1150 inertia[3] += i2 * ns_inertia_point[1];
1151 inertia[4] += -i2 * ns_inertia_point[0];
1152 inertia[5] += -i0 * ns_inertia_point[1] + i1 * ns_inertia_point[1];
1153 }
1154 } else {
1155 std::fill(inertia, inertia + 6, T(0));
1156 }
1157 }
1158
1159 if (d_inertia_d_acceleration) {
1160 std::fill(d_inertia_d_acceleration[0][0], d_inertia_d_acceleration[3][0], T(0));
1161 for (int i = 0; i != 3; ++i) { d_inertia_d_acceleration[0][i][i] = mass + ns_mass; }
1162 for (int i = 0; i != 2; ++i) { d_inertia_d_acceleration[2][i][i] = inertia_moment[i]; }
1163 d_inertia_d_acceleration[2][2][2] = inertia_moment[0] + inertia_moment[1];
1164 d_inertia_d_acceleration[2][1][0] = inertia_moment[2];
1165 if (!inertia_point_zero) {
1166 d_inertia_d_acceleration[1][0][2] =
1167 mass * inertia_point[1] + ns_mass * ns_inertia_point[1];
1168 d_inertia_d_acceleration[1][1][2] =
1169 -mass * inertia_point[0] - ns_mass * ns_inertia_point[0];
1170 d_inertia_d_acceleration[1][2][0] =
1171 -mass * inertia_point[1] - ns_mass * ns_inertia_point[1];
1172 d_inertia_d_acceleration[1][2][1] =
1173 mass * inertia_point[0] + ns_mass * ns_inertia_point[0];
1174 d_inertia_d_acceleration[2][0][0] +=
1175 mass * inertia_point[1] * inertia_point[1]
1176 + ns_mass * ns_inertia_point[1] * ns_inertia_point[1];
1177 d_inertia_d_acceleration[2][1][0] +=
1178 -mass * inertia_point[0] * inertia_point[1]
1179 - -ns_mass * ns_inertia_point[0] * ns_inertia_point[1];
1180 d_inertia_d_acceleration[2][1][1] +=
1181 mass * inertia_point[0] * inertia_point[0]
1182 + ns_mass * ns_inertia_point[0] * ns_inertia_point[0];
1183 d_inertia_d_acceleration[2][2][2] +=
1184 mass
1185 * (inertia_point[0] * inertia_point[0]
1186 + inertia_point[1] * inertia_point[1])
1187 + ns_mass
1188 * (ns_inertia_point[0] * ns_inertia_point[0]
1189 + ns_inertia_point[1] * ns_inertia_point[1]);
1190 }
1191 d_inertia_d_acceleration[2][0][1] = d_inertia_d_acceleration[2][1][0];
1192 }
1193
1194 if (d_force_and_moment_d_time) {
1195 std::fill(d_force_and_moment_d_time, d_force_and_moment_d_time + 6, T(0));
1196 }
1197
1198 if (gradient_container && compute_gradients
1199 && components == this->only_membrane_and_shear_components) {
1200 T force_moment[6];
1201 std::copy(constant_force_and_moment, constant_force_and_moment + 6, force_moment);
1202 static const int index[6] = {1, 2, 0, 4, 5, 3};
1203 for (int i = 0; i != 6; ++i) {
1204 const int ii = index[i];
1205 for (int j = 0; j != 6; ++j) {
1206 force_moment[ii] += lambda(i, j) * strain_and_curvature[j];
1207 }
1208 }
1209
1210 static const GradientContainer::FieldDescription force_beam_descr = {
1211 "SECTION_FORCES", "Fx.Fy.Fz", "Section forces", 3, 3, 1, 3, false, typeid(T)};
1212 if (gradient_container->compute_field_value(force_beam_descr.name)) {
1213 gradient_container->set_field_value(force_beam_descr, force_moment);
1214 }
1215 static const GradientContainer::FieldDescription moment_beam_descr = {
1216 "SECTION_MOMENTS", "Mx.My.dM", "Secction moments", 3, 3, 1, 3, false, typeid(T)};
1217 if (gradient_container->compute_field_value(moment_beam_descr.name)) {
1218 gradient_container->set_field_value(moment_beam_descr, force_moment + 3);
1219 }
1220
1221 static const GradientContainer::FieldDescription coor_descr = {
1222 "COOR_IP", "x.y.z", "Coordinates", 3, 3, 1, 3, false, typeid(T)};
1223 static const GradientContainer::FieldDescription disp_descr = {
1224 "DISP_IP", "dx.dy.dz", "Displacements", 3, 3, 1, 3, false, typeid(T)};
1225
1226 const T covariant_base[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
1227 MaxGradientContainer l_gradient_container;
1228 for (typename cs_properties_t::iterator p = cs_properties.begin();
1229 p != cs_properties.end(); ++p) {
1230 for (std::vector<Node*>::iterator n = p->second.begin(); n != p->second.end();
1231 ++n) {
1232 const size_t nid = (*n)->get_id();
1233 T stra[6];
1234 T stre[6];
1235 std::copy(&strain[6](0, nid), &strain[6](0, nid) + 2, stra);
1236 std::copy(&stress[6](0, nid), &stress[6](0, nid) + 2, stre);
1237 for (int i = 0; i != 6; ++i) {
1238 for (int j = 0; j != 6; ++j) {
1239 stra[j] += force_moment[i] * strain[i](j, nid);
1240 }
1241 }
1242 for (int i = 0; i != 6; ++i) {
1243 stre[0] += force_moment[i] * stress[i](0, nid);
1244 stre[3] += force_moment[i] * stress[i](1, nid);
1245 stre[5] += force_moment[i] * stress[i](2, nid);
1246 }
1247 const_cast<SolidMechanics3D<T>*>(p->first)
1248 ->get_static_linear_cross_section_gradient(
1249 covariant_base, stra, stre, &l_gradient_container);
1250 if (l_gradient_container.is_last_max()) {
1251 std::pair<int, const T*> coor = (*n)->get_coor<T>();
1252 const GradientContainer::InternalElementPosition ip = {
1253 element_coordinate, coor.second[0], coor.second[1], 1};
1254 l_gradient_container.set_current_position(ip);
1255 l_gradient_container.set_strain_and_stress(stra, stre);
1256 }
1257 }
1258 }
1259 l_gradient_container.set_gradient(gradient_container);
1260 }
1261 }
1262
1263 void field_section_integration(
1264 const T element_coordinate, const T value[6], const bool multiply_with_density,
1265 T force[6]) {
1266 if (multiply_with_density) {
1267 // Mass.
1268 force[0] = mass * (value[0] - inertia_point[1] * value[5]);
1269 force[1] = mass * (value[1] + inertia_point[0] * value[5]);
1270 force[2] =
1271 mass * (value[2] + inertia_point[1] * value[3] - inertia_point[0] * value[4]);
1272 force[3] = inertia_moment[0] * value[3] + inertia_moment[2] * value[4]
1273 + force[2] * inertia_point[1];
1274 force[4] = inertia_moment[1] * value[4] + inertia_moment[2] * value[3]
1275 - value[2] * inertia_point[0];
1276 force[5] = (inertia_moment[0] + inertia_moment[1]) * value[5]
1277 - force[0] * inertia_point[1] + force[1] * inertia_point[1];
1278
1279 // Contribution of nonstructural mass.
1280 const T i0 = ns_mass * (value[0] - ns_inertia_point[1] * value[5]);
1281 force[0] += i0;
1282 const T i1 = ns_mass * (value[1] + ns_inertia_point[0] * value[5]);
1283 force[1] += i1;
1284 const T i2 =
1285 ns_mass
1286 * (value[2] + ns_inertia_point[1] * value[3] - ns_inertia_point[0] * value[4]);
1287 force[2] += i2;
1288 force[3] += i2 * ns_inertia_point[1];
1289 force[4] += -i2 * ns_inertia_point[0];
1290 force[5] += -i0 * ns_inertia_point[1] + i1 * ns_inertia_point[1];
1291 } else {
1292 // Neutral_Axis point?
1293 UnimplementedError() << THROW;
1294 }
1295 }
1296
1297 void set_pre_stress(const T& force, const T coor[2]) {
1298 constant_force_and_moment[2] += force;
1299 if (coor) {
1300 constant_force_and_moment[0] += coor[1] * force;
1301 constant_force_and_moment[1] += coor[0] * force;
1302 }
1303 }
1304
1305 void set_pre_stress(const T force[3]) {
1306 constant_force_and_moment[2] += force[0];
1307 constant_force_and_moment[0] += force[1];
1308 constant_force_and_moment[1] += force[2];
1309 }
1310
1311 b2linalg::Matrix<T> lambda;
1312 T constant_force_and_moment[6];
1313 b2linalg::Matrix<T> v[7];
1314 b2linalg::Matrix<T> strain[7];
1315 b2linalg::Matrix<T> stress[7];
1316 Model* cs_model;
1317 typedef std::vector<std::pair<const SolidMechanics3D<T>*, std::vector<Node*> > >
1318 cs_properties_t;
1319 cs_properties_t cs_properties;
1320
1321 T neutral_point[2];
1322 T inertia_point[2];
1323 T mass;
1324 T inertia_moment[3];
1325 T ns_inertia_point[2];
1326 T ns_mass;
1327 bool compute_gradients;
1328
1329 class MaxGradientContainer : public GradientContainer {
1330 public:
1331 MaxGradientContainer() : last_is_a_maximum(false), descr_ptr(nullptr) {
1332 buffer[1] = std::numeric_limits<T>::lowest();
1333 }
1334
1335 bool is_last_max() const { return last_is_a_maximum; }
1336
1337 bool set_current_element(const Element& e) override {
1338 UnimplementedError() << THROW;
1339 return false;
1340 }
1341
1342 void set_gradient(GradientContainer* gradient) {
1343 return; // inactivated, must be fixed
1344 // if (descr_ptr == 0)
1345 // return;
1346 // gradient->set_current_position(ip);
1347 // gradient->set_field_value(*descr_ptr, buffer);
1348 // {
1349 // static const b2000::GradientContainer::FieldDescription descr = {
1350 // "SECTION_STRAINS_MATERIAL",
1351 // "Exx.Eyy.Ezz.Exy.Eyz.Exz",
1352 // (nonlinear ? "Green Lagrange strain in the material "
1353 // "reference frame" :
1354 // "Linear strain in the material reference frame"),
1355 // 3, 6, 2, 3, true, typeid(T)
1356 // };
1357 // gradient->set_field_value(descr, strain);
1358 // }
1359 // {
1360 // static const b2000::GradientContainer::FieldDescription descr = {
1361 // "SECTION_STRESSES_MATERIAL",
1362 // "Sxx.Syy.Szz.Sxy.Syz.Sxz",
1363 // (nonlinear ? "Cauchy stress in the material reference "
1364 // "frame" : "Linear stress in the material reference frame"),
1365 // 3, 6, 2, 3, true, typeid(T)
1366 // };
1367 // gradient->set_field_value(descr, stress);
1368 // }
1369 }
1370
1371 void set_current_position(
1372 const InternalElementPosition& pos, const double volume = T(0)) override {
1373 ip = pos;
1374 }
1375
1376 void set_current_volume(const double volume) override {}
1377
1378 void set_strain_and_stress(const T _strain[6], const T _stress[6]) {
1379 std::copy(_strain, _strain + 6, strain);
1380 std::copy(_stress, _stress + 6, stress);
1381 }
1382
1383 void set_field_value(const FieldDescription& descr, const T* value) {
1384 return; // inactivated, must be fixed
1385 if (value[1] > buffer[1]) {
1386 descr_ptr = &descr;
1387 std::copy(value, value + descr.nb_comp, buffer);
1388 last_is_a_maximum = true;
1389 } else {
1390 last_is_a_maximum = false;
1391 }
1392 }
1393
1394 bool compute_field_value(const std::string& name) const override { return true; }
1395
1396 private:
1397 bool last_is_a_maximum;
1398 const FieldDescription* descr_ptr;
1399 InternalElementPosition ip;
1400 T buffer[20];
1401 T strain[6];
1402 T stress[6];
1403 };
1404};
1405
1406} // namespace b2000
1407
1408#endif // B2REISSNER_BEAM_PROPERTY_H_
THROW
#define THROW
Definition b2exception.H:198
b2000::GradientContainer
Definition b2solution.H:54
b2000::GradientContainer::set_field_value
virtual void set_field_value(const FieldDescription &descr, const int *value)
Definition b2solution.H:146
b2000::RTable
Definition b2rtable.H:427
b2000::ReissnerBeamElasticCProperty
Definition b2reissner_beam_property.H:921
b2000::ReissnerBeamElasticIProperty
Definition b2reissner_beam_property.H:866
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314